Calhoun: The NPS Institutional Archive

Theses and Dissertations Thesis Collection

1985

The fractal geometry of nature; its mathematical

basis and application to computer graphics.

Gaddis, Michael E.

http://hdl.handle.net/10945/21349
DUDLEY KNOX LIBRARY
NAVAL POSTGRADUATE SCHOOL
MONTEREY, CALIFORNIA 93943-8003
NAVAL POSTGRADUATE SCHOOL
Monterey, California

THESIS
The Fractal Geometry of Nature;
Its Mathematical Basis and
Application to Computer Graphics
by
Michael Edward (Caddis
December 19S5

Thesis Advisor: M. .1. Zvda

Approved for public release: distribution i- unlimited.

T226311
CURITY CLASSIFICATION OF THIS PAGE

REPORT DOCUMENTATION PAGE

REPORT SECURITY CLASSIFICATION lb. RESTRICTIVE MARKINGS
UNCLASSIFIED
SECURITY CLASSIFICATION AUTHORITY 3 DISTRIBUTION/AVAILABILITY OF REPORT
Approved for public release;
DECLASSIFICATION / DOWNGRADING SCHEDULE distribution is unlimited
PERFORMING ORGANIZATION REPORT NUMBER(S) 5. MONITORING ORGANIZATION REPORT NUMBEP(S)

NAME OF PERFORMING ORGANIZATION 6b OFFICE SYMBOL 7a. NAME OF MONITORING ORGANIZATION

(If applicable)
aval Postgraduate School Naval Postgraduate School
ADDRESS (City, State, and ZIP Code) 7b. ADDRESS (City, State, and ZIP Code)

onterey, CA 93943-5100 Monterey, CA 93943-5100

NAME OF FUNDING /SPONSORING 8b. OFFICE SYMBOL 9 PROCUREMENT INSTRUMENT IDENTIFICATION NUMBER
ORGANIZATION (If applicable)

ADDRESS (City, State, and ZIP Code) 10 SOURCE OF FUNDING NUMBERS

PROGRAM PROJECT TASK WORK UNIT
ELEMENT NO NO NO ACCESSION NO

TITLE (Include Security Classification)

HE FRACTAL GEOMETRY OF NATURE; ITS MATHEMATICAL BASIS AND APPLICATION' TO

OMPUTER GRAPHICS
PERSONAL AUTHOR(S)
ichael E. Gaddis
a TYPE OF REPORT 13b TIME COVERED 14 DATE OF REPORT (Year, Month, Day) 15 PAGE COUNT
Master's Thesis FROM TO 1985 December 151
SUPPLEMENTARY NOTATION

COSATI CODES 18 SUBJECT TERMS (Continue on reverie if necessary and identify by block number)
FIELD GROUP SUB-GROUP Computer Graphics, Fractals, Stochastic Terrain
Modelling, Terrain Modelling

ABSTRACT (Continue on reverse if necessary and identify by block number)

ractal Geometry is a recent synthesis of old mathematical constructs. It
as first popularized by complex renderings of terrain on a computer graph
iedium. Fractal geometry has since spawned research in many diverse scien ti
"ic disciplines. Its rapid acceptance has been achieved due to its abilit y
o model phenomena that defy discrete computation due to roughness and dis
.ontinuit ies With its quick acceptance has come problems.
. Fractal ;^eome I

s a misunderstood idea that is quickly becoming buried under grandiose

:erminology that serves no purpose. Its essence is induction using simple
;eometric constructs to transform initiating objects. The fractal objects
hat we create with this process often resemble natural phenomenon.
iurpose of this work is to present fractal geometry to the graphics pro m
ler as a simple workable technique. We hope to demystify the concepts of
fractal geometry and make it available to all who are interested.
o distribution /availability of abstract 21 ABSTRACT SECURITY CLASSIFICATION
Ol'MCLASSIFIED/UNLIMITED SAME AS RPT D DTIC USERS UNCLASSIFIED
2aNAMF OF RESPONSIBLE INDIVIDUAL 22b TELEPHONE (Include Area Code) 22c
'rof. Michael J. Zyda 408-646-2305 5 2 Zk

DFORM 1473, B4 mar 33 APR edition may be used until exhausted SECURI t y CLASSIFICATION OF ~~»s ?ACi
Ail other editions are obsolete
1
 Approved for public release, distribution unlimited

The Fractal Geometry of Nature;

ItsMathematical Basis and
Application to Computer Graphics
by

Michael Edward £Jaddis

Captain, United States Marine Corps
B.S.B.A, University of Florida, 1978

Submitted in partial fulfillment of the

requirements for the degree of

MASTER OF SCIENCE IN COMPUTER SCIENCE

from the

NAVAL POSTGRADUATE SCHOOL

December 1985
 A BSTRAC T
Fractal Geometry is a recent synthesis of old mathematical constructs. It was
first popularized by complex renderings of terrain on a computer graphics
medium. Fractal geometry has since spawned research in many diverse scientific

disciplines. Its rapid acceptance has been achieved due to its ability to model
phenomena that defy discrete computation due to roughness and discontinuities.

With its quick acceptance has come problems. Fractal geometry is a

misunderstood idea that is quickly becoming buried under grandiose terminology

that serves no purpose. Its essence is induction using simple geometric constructs

to transform initiating objects. The fractal objects that we create with this

process often resemble natural phenomenon. The purpose of this work is to

present fractal geometry to the graphics programmer as a simple workable

technique. We hope to demystify the concepts of fractal geometry and make it

available to all who are interested.

3
 5
1
TABLE OF CONTENTS

I. AN INTRODUCTION TO FRACTAL GEOMETRY 7

A. MATHEMATICS AS A MODEL FOR OUR UNIVERSE 7

B. FRACTAL GEOMETRY 9

C. GOALS OF THIS RESEARCH 9

II. THE MATHEMATICAL BASIS OF FRACTAL

GEOMETRY 10

A. PRELIMINARIES 10

B. DIMENSION 15

C. FRACTAL CURVES AND SETS 24

III. THE IMPLEMENTATION OF FRACTALS IN COMPUTER

GRAPHICS 35

A. THE IMPLEMENTATION PROBLEM 35

B. MAPPING FRACTALS TO A BOUNDED SPACE 38

IV. FRACTAL COMPUTATION IN R 2

45

A. THE GEOMETRY OF INITIATOR - GENERATOR 45

B. THE MID-POINT DISPLACEMENT TECHNIQUE 49

C. A KOCH-LIKE FRACTAL ALGORITHM 51

D. IMPLEMENTATION STRATEGIES 55

E. SUMMARY 55

V. FRACTAL GEOMETRY FOR GRAPHICS TERRAIN 61

A. MODELING MOUNTAINOUS TERRAIN 61

B. FRACTAL TOOLS FOR TERRAIN MODELING 68

VI. SHORT CUTS TO MOUNTAIN SHAPES 87

A. RECTANGULAR MIDPOINT TECHNIQUE 87

B. PARAMETRIC CUBIC SURFACES 02

VII. CONCLUSIONS 100

A. DIRECTIONS FOR FURTHER STUDY 100

B. CONCLUSIONS 102

APPENDIX A: FRACTAL COMPUTATION IN R 2

103

APPENDIX B: RANDOM NUMBER GENERATORS 112

APPENDIX C: THE TRIANGULAR MOUNTAIN 116

APPENDIX D: THE RECTANGULAR MOUNTAIN 121

APPENDIX E: GEOMETRIC SUPPORT 126

LIST OF REFERENCES 128

INITIAL DISTRIBUTION LIST 129

 ACKNOWLEDGEMENTS
I wish to express my gratitude to a number of people who supported this

study. To my advisor, Dr. Michael Zyda, who gave me the opportunity to

pursue a subject that was both interesting and challenging. Who also allowed me
considerable leeway in my research and had the good sense to keep me on track

when I strayed too far. To Leut. Steven Firth, RAN, for Ozdraw, his

excellent figure generation system, which made the creation of figures for my
thesis a joy. To Lt. John Ynrchak, USN, who was the intrepid hacker that

developed the thesis macros for the QMS laser printer. But most of all I would
like to thank my wife, Kathleen, who has shown infinite patience with my
seemingly endless late hours and far-off gazes. Who has given me a wonderful

daughter during this time, and who gives me the joy of life.

The measure of any man is his family.

6
 I. AN INTRODUCTION TO FRACTAL GEOMETRY

A. MATHEMATICS AS A MODEL FOR OUR UNIVERSE

The various branches of mathematics have through time developed as a

response to the need for more detailed models to describe new developments,

both technological and philosophical. This was true when Newton developed

calculus and also true during the late 1800's through the 1920's when a schism

developed between the classical mathematicians and some brilliant innovative

thinkers.

1. The Mathematical Crises of the Early 19th Century

One of man's greatest strengths is his ability to question his surroundings

and beliefs and through this questioning develop new insight and innovation.

Most mathematical systems are developed for use in applications. Man's natural
inquisitiveness often leads him to develop his systems beyond the application and

into abstract theory. This theory drives him to investigate the applications and

often yields direction for new discoveries that were not previously foreseen or that

defy intuition.

Georg Cantor (1845-1918) was the most notable of a number of

mathematicians who questioned the basic precepts of mathematics and developed

the modern set theory. Some of Cantor's discoveries seemed to invalidate many
of the long held beliefs of mathematics. Cantor and his peers became deeply
involved in controversy over their findings. Their discovery of functions which

seemed to violate the basic rules of geometry and calculus were deemed as

monsters and unworthy of consideration by reasonable men because they lacked

usefulness to any application then known [Ref. l:pp. 9]. These new concepts

would remain in the arena of pure theoretical mathematics until science

developed to a point where the old models could no longer adequately describe it>

processes and would look to the new mathematics for a new perspective.

It was from these discoveries that Fractal Geometry was born Ref.

l:Chap. 2]. It will be seen in the following chapters that fractal geometry
is a synthesis of many of the concepts which developed from the mathematical

schism of the 19th century, most notably set theory and topology.

2. What is a Mathematical Model

Reference 2 defines a mathematical model in the following fashion [Ref.

2:pp 1-3]:

A mathematical model is a mathematical characterization of a phenomenon or

process. It has three essential parts: a process or phenomenon which is to be
modeled, a mathematical structure capable of expressing the important proper-
ties of the object to be modeled, and an explicit correspondence between the
two.

Although the phenomenon of interest need not be taken from the real world,
they usually are. The real world component is described quantitatively by such
things as parameter values and at which time things occur.

The second component of a model is an abstract mathematical structure. In it-

self, the structureabstract and has no intrinsic relation to the real world.
is
However because of its abstractness it can be used to model many different
phenomena. Every mathematical structure has an associated language for mak-
ing assertions. If the mathematical model is successful, the language of its
mathematical structure can be used to make assertions about the object being
modeled.

The third component of a model is a specification of the way in which the real
world is represented by the mathematical structure, that is, a correspondence
between the elements of the first component and those of the second.

3. The Euclidean Model

When using mathematics to describe man-made objects, the Euclidean

model (standard Euclidean geometry) is usually satisfactory. Its structure is

simple and pure, which appeals to an engineer's nature. But as technology

expands and we need to describe processes that are not well behaved, we need to

develop a geometry that can adequately model our process within a certain

closeness of scale.

No model can completely describe a natural object because nature does

not follow the man-made rules that we impose on our model. But at a given

scale, the model (if it is accurate) can describe the object with enough precision

to be of help in constructing it. Engineers use the geometry of a straight line to

describe a wall but this wall, when viewed closely enough, is not straight at all.

This is of no matter to the engineer because his model is accurate for his scale of

reference.

8
B. FRACTAL GEOMETRY
One man who saw a need for a new geometry was Benoit B. Mandlebrot. He
felt that Euclidean geometry was not satisfactory as a model for natural objects.

To anyone who has tried to draw a picture of a nonregular object (such as a tree)

on a computer graphics screen, using the Euclidean drawing primitives usually

provided, this is an obvious statement. The strength of Mandlebrot's finding was
his research into the findings of the earlier mathematicians and the development
of a practical application of their theory. Mandlebrot coined the term Fractal to

describe a class of functions first discovered by Cantor (Cantor's dustj. Koch (the

Koch curve) and Peano. He showed how these functions yield valuable insight

into the creation of models for natural objects such as coastlines and mountains.
Mandlebrot popularized the notion of a fractal geometry for these types of

objects. Although he did not invent the ideas he presents, Mandlebrot must be

considered important because of his synthesis of the theory at a time when

science was reaching out for new more accurate models to describe its processes.

C. GOALS OF THIS RESEARCH

There are two approaches that can be taken in the investigation of fractal

geometry and computer graphics.

- To view the computer as a tool to enhance the investigation of fractal

geometry.
or
- To view fractal geometry as a tool to enhance the realism of computer
graphics.

1
This research will take the later approach . It is designed to investigate the

mathematics of fractal geometry and to show its application to computer

graphics. I hope to be able to tame the subject of fractal geometry by making its

mathematics and technique accessible to the average computer scientist.

1
Where Mandelbrot took the former.

9
 II. THE MATHEMATICAL BASIS OF FRACTAL GEOMETRY

This chapter is a brief introduction to the mathematical foundations that

underlie the theory of fractals. Little technique currently exists for the practical

application to attain complete mathematical rigor when using fractal functions

(i.e. it is very difficult to prove that a set is fractal). This causes the non-

mathematician to accept much of what he does with fractals on faith. It is

instead important to understand the theory intuitively. This can be gained by a

cursory look at the mathematical foundations for fractals.

A. PRELIMINARIES
A complete definition of fractals is given later in this chapter but before we
can understand that definition, we must establish a foundation in set theory.

Fractals were discovered in set theory and topology. They can be considered as

an outgrowth of investigations into these related fields.

1. What is a Set

A set is defined in [Ref. 3:pp. 11] in the following fashion:

A set is formed by the grouping together of single objects into a whole. A set is
a plurality thought of as a unit. We can consider these statements as expository,
as references to a primitive concept, familiar to us all, whose resolution into
more fundamental concepts would perhaps be neither competent nor necessary.
We will content ourselves with this construction and will assume as an axiom
that an object M
inherently determines certain other objects a, b, c, ... in some
undefined way, and vice versa. We express this relation by the words: The set
consists of the objects a, b, c, ...
M
This definition is intentionally vague to allow the set to become the basic

building block for all mathematical constructs.

2. Some Set Theoretic Concepts

This section presents some background definitions for concepts used in

the body of this chapter. The reader is directed to the references for a detailed

explanation or proof [Ref. 2:Chap 5].

Definitions:

Cardinality
Two sets S and T are said to have the same number of elements, or to have the
same cardinality, if there is a one-one function /from 5 to T.

10
 Finite and Infinite Sets

A Set 5 issaid to be finite if 5 has the same cardinality as 0, or if there is a po-

sitive integer n such that S has the same cardinality as {1,2,3,4, 5, ... n }. Other-
wise S is said to be infinite.

Countability
A set S is said to be countable if S has the same cardinality as a subset of N, the
set of positive integers. Otherwise, 5 is said to be uncountable.

Propositions:

a). Any subset of a finite set is finite.

b). Any subset of any countable set S is countable.

c). The set of natural numbers N is countable.

d). The set of rational numbers Q is countable.

d). The set of irrational numbers is countable.

e). The union of a countable collection of countable sets is countable.

f). The set of real numbers R is uncountable.

3. Some Topological Concepts

This section presents some topological background concepts used in the

body of this chapter. The reader is directed to the references for a detailed

explanation [Ref. 3].

Concepts:

Metric Space
A metric space is a set in which we have a measure of the closeness or proximity
of two elements of the set, that is, we have a distance defined on the set. For ex-
ample, a metric on R 2
would be the pythagorean metric:
2 2
D((xi,yi),(x 2 ,y 2 )) = x/(x 1
-i 2 ) +(yi-J/2)

Covering a Set
Let A' be a topological space and S a subset of X. A cover of the set 5 is exactly
what its name implies, a collection of subsets of X which cover S, that is, whose
union contains S.

4. What is a Function

A function is defined in [Ref. 2:pp. 193-194] as

A function from a set A to a set B is a rule which specifies an element of B for

each element of A.

Definition:

Let A and B be sets. A function (or map, or transformation) f from A to B.

denoted / A -» B is a relation from A to B such that for every a which is an
: .

element of A, there exists a unique b in B such that < a,b> is an element of /.

We write f(a) = b.

11
 If /is function from A to B, then A is called the domain of the function /and B
is called the codomain of /.

To completely define a function we must specify the domain, the codomain, and
the value f(x) for each possible argument x.

Functions can be viewed as a specification of a method to describe the

creation of a set from other sets using some agreed upon mathematical
symbolism. The functions can yield powerful results when the target set (co-

domain) is complex and not easily described by set theoretical constructs. This is

especially true in fractal functions when the domain is R 1

and the object created
(set, co-domain) is a nonregular shape. This is one reason why the computer

graphics system is useful in the investigation of fractal functions. The computer

can model the infinite function and display a finite approximation of the created

fractal set.

5. Useful Functional Concepts

It is often helpful to clearly understand the universe of discourse within

which a function exists. The function can be rigorously defined within the above

constructs but lack intuitive appeal due to its complexity. Mathematicians have

defined many useful concepts to describe functions. The concepts applicable to

this study are described below.

a. Partial Functions

Most of the fractal functions in this study have as their domain some
undefined subset of R . It is useful then to consider them as partial functions

and not concern ourselves with a rigorous description of the domain of the

fractal. We take our definition from [Ref. 2:pp. 201-202].

It is often convenient to consider a function from a subset A ' of A to a set B

without exactly specifying the domain A of the function. Alternatively, we can'

view such a situation as one where a function has a domain A and codomain B,
but the value of the function does not exist for some arguments of A. This is
called a partial function.

Definition:

Let A and B be sets. A partial function f with domain A and codomain B is any
function from A' to B, wnere A' is a subset of A. For any x which is an element
of A - A', the value of / (x is said to be undef ined ) .

b. Bijectivity

It is often useful to know to what extent a function maps from the

domain to the codomain. If a function is not a partial function and every point of

12
the domain maps to a point in the codomain then we want to know if all points

of the domain A in the mapping /(A) map to distinct points in the codomain B.
We may also want to know to what extent the mapping /(A) covers the set B.

The definition of bijective, surjective and injective functions is from [Ref 2:pp.
204].

Definition:

Let / be a function /: A »B.

(a) / is surjective (onto) if /(A) = B,

(b) / is injective (one-to-one) if a / a '

implies /(a) # /(a'),

(c) / is bijective (one-to-one and onto) if / is both surjective and injective.

6. Functions From RN - RN
A point in R' space is specified by an n- tuple of the form

(jr
1
,x 2
,x 3 ,
xn ). To completely specify a function from R^ — R N each point in

the domain must map to a point in the codomain. An example is:

/ : R -R 2 2

2
/((x 1? x 2 ))= (x, ,,/)

This function is well defined. For each point in the domain of the function we
have specified a unique point in the codomain.

Most of the functions that are covered in this study are mappings within
R 2
or R3 . Fractal sets exist in all finite dimensions but it is impractical at this

point to use fractal functions beyond the fourth dimension in view of the graphics

display medium's limitation to two dimensions. The use of fractal functions

whose dimension is between 3 and 4 is currently being investigated by allowing

the function to roam the fourth dimension and then taking time slices which

yield three dimension approximations of the set [Ref. 4].

7. Inductive Definitions of Sets

Functional constructs do not always provide a convenient means of

charactering an infinite set. It is sometimes more eloquent and powerful to use

the inductive method to characterize a set.

13
 Our definition of inductive definitions of sets is from [Ref. 2:pp. 199-201].

An inductive definition of a set always consists of three distinct components.

1. The basis, or basic clause, of the definition establishes that certain objects
are in the set. The basic clause establishes that a set is not empty and charac-
terizes the "building blocks" (the seeds of the induction) which are used to con-
struct the set from the inductive clause.

2.The induction, or inductive an inductive definition establishes

clause, of
the ways in which elements of the can be combined to obtain new elements.
set
The inductive clause always asserts that if objects z ,y ,...,z are elements of the
set. then they can be combined in certain specified ways to create other elements
of the set (thus from the basic clause (or seeds) of the induction we induce the
remaining elements of the set).

3. The extremal clause asserts that unless an object can be shown to be a

member of a set by applying the basis and inductive clauses a finite number of
times, then the object is not a member of the set.

An example of an inductively defined set is:

(Basis)

t A
(Induction)

If n e A, then (n +2) e A
(Extremal)

No integer is an element of A unless it can be shown to be so in a finite number

of applications of clauses 1 and 2 above.

The set that we defined is the set of all even nonnegative integers.

8. The Path To Fractals

The path to fractals by the non-mathematician is not through theory but

through the investigation of their functions and methods of construction. This

investigation (and experimentation) yields considerable insight into the nature of

fractals.

The choice of which set-descriptive methodology to use (functional or

inductive) in describing a set is often a matter of style but can be dictated by

necessity if one method is inordinately tedious.

Most of the fractal functions that are introduced in this study use the

inductive method as the primary functional tool. In fact, these functions are a

hybrid of the functional and inductive constructs described above.

Often called a recurrence definition.

14
B. DIMENSION
The classification of fractal sets from non-fractal sets is based on the

dimensional qualities of the set. To understand fractals you must have an

appreciation for these differences.

The concept of dimension is one rife with difficulties. Many of the gr^at

mathematicians have attempted to define dimension as a rigorous concept

consistent with the known mathematical systems. For each of their attempts

however, the concept becomes more prone to contradiction and paradoxes.

There currently exist five definitions of dimension that date back to the late

1800V. The classification of Fractal sets into a class of sets is the result of the

discovery of functions that created sets which did not fit comfortably into the

topological definition of dimension (which was the accepted definition at the

time). Fractal sets are rigorously classified as those sets that demonstrate a

difference between the Hausdorff-Besicovitch dimension and the standard

topological dimension.

1. An Intuitive Approach to Dimension

Dimension is a concept that seems intuitive when it is first introduced in

Euclidean geometry as the standard three dimensions. Long after Euclid made
the first attempts at defining dimension and concurrent with the discovery of

atomic particle physics the concept of dimension was rethought by the prominent

mathematicians of the time. This was necessary to realign the mathematical

model for the geometry of objects with the new view of what those objects were

made of. Our increasing ability to focus on the nature of matter inevitably

causes the models we use to change.

The dilemma that arose from the new concepts of dimension quickly

developed into a theoretical debate that left intuition behind. When human
intuition fails, we must rely upon well founded models that are based on axioms

of basic mathematical truth. It is only through the rigorous investigation of our

mathematical models that allows us to go beyond intuition and investigate the

1
1). Cantor and Minkowski; 2). Bouligand and Minkowski; 3). Pontrajgin. Schnirelman and
Kolornogorov, Tihomirov; 4). Hausdorff-Besicovitch and 5). the topological dimension (there are
others). Most of these definitions are concerned with the most efficient method of covering a set
(i.e. are topological concerns).

15
true physical nature of the objects we model. The debate still rages today and

borders on the philosophical. Two examples should suffice to demonstrate the

complexity and possible paradoxes that can arise from dimension theory.

The first is from [Ref. 5:pp. 323-344].

Consider the way in which we define the density of air at a given point and at a
given moment. We picture a sphere of volume V
centered at that point and in-
cluding the mass Wi. The quotient M/V is the mean density within the sphere,
and by the true density we denote some limiting value of this quotient. This no-
tion, however, implies that at the given moment the mean density is practically
constant for spheres below a certain volume. This mean density may be notably
different for spheres containing 1,000 cubic meters and 1 cubic centimeter
respectively.

Suppose the volume becomes continually smaller. Instead of becoming less and
less important, these fluctuations come to increase. For scales at which the
brownian motion shows great activity, fluctuations may attain 1 part in 1,000,
and they become of the order of 1 part in 5 when the radius of the hypothetical
spherule becomes of the order of a hundredth of a micron.

One step further and our spherule becomes of the order of a molecule radius. In
a gas it will generally lie in intermolecular space, where its mean density will
henceforth vanish. At our point the true density will also vanish. But about
once in a thousand times that point will lie within a molecule, and the mean
density will be a thousand times higher than the value we usually take to be
the true density of the gas.

Let our spherule grow steadily smaller. Soon, except under exceptional cir-
cumstances, it will Decome empty and remain so henceforth owing to the intra-
atomic space; the true density vanishes almost everywhere, except at an infinite
number of isolated points, where it reaches an infinite value.

The second is from [Ref. l:pp. 17-18].

Consider a ball of 10 cm diameter made of a thick thread of 1 mm diameter that

(in latent fashion) possesses several distinct effective dimensions.

To an observer placed far away, the ball appears as a zero-dimensional figure: a

point. of 10 cm resolution, the balj of thread is a three-
As seen from a distance
dimensional figure. At 10 mm, it is a mess of one-dimensional threads. At 0.1
mm, each thread becomes a column and the whole becomes a three-dimensional
figure again. At 0.01 mm, each column dissolves into fibers, and the ball again
becomes one-dimensional, and so on, with the dimension crossing over repeated-
ly from one value to another. When the ball is represented by a finite number of
atomlike pinpoints, it becomes zero-dimensional again.

It is interesting to note that each of these examples demonstrate

dimension as a reflection of physical properties dependent on the observers point

of reference. Each ends with reference to the paradox of atomic particles. That
paradox is, for any collection of finite (or countably infinite) points, the

dimension is zero [Ref. 6:pp. 1-8]. Since the earth and sun each have a finite

collection of atoms then accordingly their dimension is zero. Dimension exists

16
only for a mathematical continuum and as such lacks application to the physical
4
universe as we currently know it .

The two properties (continuity and dimension) cannot be separated.

Before the advent of atomic theory, matter was viewed as continuous and

composed of basic elements that were indivisible. While the debate raged over the

practical and philosophical aspects of the nature of matter it became apparent

that the mathematical models which represent matter would have to change. It is

not practical to represent objects by representing each atom and its position

relative to the entire set. The power of modeling would thus be lost; that is. the

ability to model complex objects and their interaction by relatively simple

constructs. Thus the fact is reinforced that models can only represent objects

through gross approximations and that the model is only effective for a restricted

frame of reference. Without this realization, dimension would have very little

application.

2. Topological Dimension

a. An Intuitive Approach
The concept of dimension is very old. It is based on the algebraic

concepts of Euclidean n space and the notion that a set has dimension n if the

least number of real parameters needed to describe its points was n. This fuzzy

definition was accepted for a very long time until the advent of Cantor and set

theory. Cantor showed that dimension can be changed by a 1-1 transformation

from an interval to a planar object. The fuzzy notion of dimension, as defined,

was challenged and required rethinking.

The mathematicians who did not accept many of the findings of set

theory at the time (but who could not disregard Cantors findings) began to

consider ways of explicitly defining dimension. The new definition would have to

be applicable to the bizarre functions of Cantor. Koch and Peano as well as the

4
The set theoretical concepts of finiteness, countably infinite and uncountabl) infinite (con-
tinuous) sets carry with them very profound implications. It is premature to view matter a> mere-
ly collections of finite atoms. Science may yet find true continuity (in the mathematical sense) in
atomic matter and the universe. For now, matter is what it is and our pronunciations upon it will
not change its true texture.

17
relatively simple objects that had previously fit into the old definition without

contradiction.

There was one crucial problem that Cantor's findings raised

[Ref. 6:pp. 4-6]:

An extremely important question was left open: Is it possible to establish a

correspondence between Euclidean n space ana Euclidean m space combining the
features of both Cantor's and Peano's constructions, i.e. a correspondence which
is both 1:1 and continuous? The question is crucial since the existence of a
transformation of the stated type between Euclidean n space and Euclidean m
space would signify that dimension (in the natural sense that Euclidean n space
has dimension n ) has no topological meaning whatsoever.

This fundamental problem was answered in 1911 by Brouwer. He

proved that Euclidean n space and Euclidean m space were not homeomorphic

unless n equals m. To say that two spaces A and B are homeomorphic means
that a mapping / :A^B exists, such that /is continuous over .4 and bijective.
~
Additionally, the inverse of this mapping /
l
B^A , is continuous over B and
bijective. If two spaces are homeomorphic then it is analogous to saying that they

are topologically equivalent.

Further research was done and a precise definition of topological

dimension of a set was developed. This definition assigned an integer value as the

dimension of any set based on its topological properties,

b. Definition of Topological Dimension n

The rigorous investigation of dimension is beyond the scope of this

study but the following definition is included for completeness [Ref. 6:pp. 24]:

Roughly speaking, we may say that a space has dimension ^ n if an arbitrarily

small piece of the space surrounding each point may be delimited by subsets of
dimension < n - 1. This method of definition is inductive, and an elegant start-
ing point for the induction is given by prescribing the null set as the (-1) dimen-
sional space.

Definition:

The empty set and only the empty set has dimension -1.

A space A' has dimension ^ n ( n > 0) at a point p if p has arbitrarily small

neighborhoods whose boundaries have dimension ^ n - 1.

X has dimension ^ n, dim X ^ n, if .V has dimension < n at each of its points.

X has dimension n at point p if it is true that X has dimension -^ n at p and it is

false that X has dimension ^ n - 1 at p.
X has dimension n if dim X< n is true and dim X ^ n - 1 is false.

X has dimension oo if dim X^ n is false for each n.

18
 The topological dimension is rigorous and consistent for all sets that

exist within a metric space. The problem that arises with fractal sets and its

topological dimension is not that the topological dimension is wrong. Fractal sets

like all sets in a metric space exhibit a topological dimension. The question is

then, is the topological dimension an accurate description of the dimension of the

set or can we find a better way to characterize the dimensional qualities of the

set? This question can be extended; is the topological definition of dimension

useful and consistent with the notion of dimension and space? Can we devise a
3
better concept which can further refine dimension and make it more useful .

3. The Hausdorff-Besicovitch Dimension

This section is intentionally brief due to the subject's complexity and to

the lack of practical technique that it yields. The Hausdorff measure of a set is a

complex characterization of a method for covering a set. Hausdorffs theorem is

proved using the existential qualities of infinite sets in a metric space H. While

the theorem may be important to mathematical theory, it proves unfortunate

that there is no straightforward practical method for determining the Hausdorff

measure of a set.

a. An Intuitive Approach to the Hausdorff Dimension

The acceptance of the Hausdorff method for covering a set as a

measure of dimension is not universal [Ref. 6:pp. 102] and [Ref. l:pp. 363-365].

The debate is between the disciplines of topology and metrics and is not wholly

germane to this study. It is beneficial to divorce ourselves from the debate and

consider both the topological dimension DT and the Hausdorff-Besicovitch

dimension (HB) as merely measures of qualities of a set's structure. Certainly

sets exist that have a topological dimension equal to 1 but in no way resemble a

simple Euclidean curve If the Hausdorff dimension yields a better measure of a

set's structure that provides a mathematical and intuitive difference that is useful

to us, it would be beneficial for us to investigate it.

5
Try not to ascribe grandiose implications to a set's dimension (the fourth dimension as time
or some such) as this is premature at best. Rather, view a set's dimension as merely descriptive
terminology the terminology of bijectivity describes a function's characteristics The
much like
problem most people have with this mental abstraction is the visual reinforcement that the\ re-
ceived from the notion of the standard three dimensions.

19
 The Hausdorff measure of a set was developed during the same
period that the new topological dimension was invented to solve the paradoxes of

Cantor. The topological dimension was based on the idea of a neighborhood of a

point within a Euclidean space of R . The connection to metric spaces and the
idea of measure is obvious when you consider that the Hausdorff measure of a set

is also based on this notion of a spherical neighborhood and what Hausdorff calls

the test function of a set. The test function of a set denoted h(p ) is a function

that characterizes the "best" method of covering a point with the spherical ball

of radius p that covers points of the set, which in their union, cover the entire

set.

Consider for example the test function for a surface within R\ A

surface can be covered by discs (circles). The formula for the area of a circle

becomes the test function for the surface. The formula for the area of a circle

always contains the constant factor it multiplied by the square of the radius r.

This radius is the measure p as above. This leaves us with a test function for a

planar shape in R 2
of h(p )
= wp 2 .

You might expect that the test function for a spherical neighborhood

in a Euclidean space above R3 would be very difficult to imagine, and indeed it

is. Hausdorff further complicated the idea of test functions (even within the lower

dimensions) by allowing a test function to assume a non-integer parameter d so

that the test function h(p )

= 0(d)p could have a real-valued parameter d.

This further refinement of the test function allowed Hausdorff to make assertions

about how this test function h(p ) behaved when the parameters p and d were

allowed to vary.

Hausdorff imagined the parameter p reducing in size until it

approached zero. The effect of this on our disc example is increasingly smaller

and smaller discs around points of the planar set. As the disc size is decreased,

fewer points of the set are contained in each disc neighborhood. This requires

more discs to cover the set. As the parameter p becomes infinitely small the

number of discs required to cover the planar set approaches oc. We allow the

parameter to grow arbitrarily small. It is interesting to study the test function

20
and see what happens to the total area when the areas of the collection of discs

which cover the planar set are summed.

Let's reflect upon the mathematical process that we are developing.

When we attempt to approximate the area of the planar set by the union of the
discs which make up its cover, we are essentially observing small patches of the

surface and approximating the area of the set by making assertions about the
intrinsic qualities of the patch. The notion that this measure is merely an

approximation is important. As the size of our patch grows increasingly small, we

can expect that we will get a better fit with our patches and hence a better

approximation of the area. The notions of approximation and fit become

increasingly helpful when you consider functions which describe sets of infinitely

rough texture as we find in fractal sets.

The importance of Hausdorffs discovery lies in the fact that for a

test function of a set, the parameter d is special. As p-»0 he discovered that

there existed a unique real number d such that for d'<d. the infinum defined

by the test function using d' approaches oc (for any countably infinite set). And
for a d'>d the corresponding infinum approaches zero. This means that every

set has a parameter which can be associated with it that is closely related to the

amount of space that it occupies. This number is the Hausdorff-Besicovitch

dimension (measure) of the set.

These results only tell us that a number and function exist. They do

not tell us how to compute them in the general case. This gives us the quandary

of dealing with the HB dimension as a known concept that we can make

allusions to, but can rarely compute (at least at the present time),

b. Definition of the Hausdorff Dimension

The formal definition of the Hausdorff dimension requires that we

formalize the inttiitive discussion above. We first define what a ;>-measure of a set

is (analogous to the disc above) [Ref. 6:pp. 102-103].

A p-dimensional measure for each non-negative real number p was defined bv

Hausdorff for arbitrary metric spaces. This measure is a metrical concept. whiTe
dimension is purely "topological. Nevertheless there is a strong connection
between the two concepts, for it turns out that a space of [topological] dimen-
sion n must have positive rt-dimensional measure (the Hausdorff measure).

For our example this would be the sum of the areas of the discs

21
 Definition:

Let X be a space and p an arbitrary real number. ^ p <oo. Given e >0 let

m< p = inffj [6{A,)\>

1 = 1

where X = A + A + A* +
1 2
• •
is any decomposition of X into a countable
number of subsets of diameter less than e , and the superscript p denotes ex-
ponentiation. Let

m. (X) = sup m% (X).

t >o

mp (X) is called the p-measure or (p-dimensional) measure of X.

Proposition:

If p <q then m p
(X) ^ m ?
(X); in fact p <q and mp (X) <oo imply m ?
(X) = 0.

Conversely, if p >q then mp (X) ^ m q

(X); and if p >g and
mp (X) <oc implies m q
(X) = oo .

c. Mandelbrot's Misgivings about the HB Measure

It is clear from the previous definition that for practical applications

the Hausdorff dimension is difficult to compute directly. In [Ref. l:pp. 14-19] and

[Ref 7], Mandelbrot expressed a distaste for the focusing of attention on the HB
dimension. He states:

''I developed the definition of fractals using the topological and Hausdorff di-
mensions in response to colleagues who urged me to do so. They felt that it was
necessary to rigorously define the concept within firm mathematical criteria. I
have come to believe that an empirical definition would be more beneficial at
this time because the present definition denies the inclusion of some shapes that
could best be described as fractals."

Mandelbrot believes that the definition of the Hausdorff dimension is

too difficult to deal with and is perhaps too restrictive. He prefers to focus on the

behavioral aspects of the recurrence relationship involved in fractals and the

empirical results from these equations.

A practical application of fractal functions in computer graphics

does, by necessity, bend to this same paradigm. This realization should not blind

us to the fundamental nature of a fractal equation's uniqueness, however. It is

important to understand the dimensional aspect of the fractal discourse to

appreciate the potential importance of fractals.

It is not preordained that fractal equations model nature with a

greater degree of accuracy then does the Euclidean model. The future may prove

the fractal model the superior method, however. The dimensional qualities of

fractal functions may be the aspect that proves this to be so.

22
 4. Why Consider Dimension

For the purpose of this study, it is not important that a rigorous feel for

the mathematical properties of dimension theory be grasped. In fact one needs no

knowledge of dimension to use fractal techniques in the generation of computer

graphics terrain. The literature is rife with articles about fractal objects in

computer graphics and it seems de rigueur to include an approximate fractal

dimension as part of its description. The techniques used to approximate

dimension as presented in [Ref. l:pp. 56-57] are mathematically unproven'. More

importantly, the dimension yields little intuitive insight: one is hard pressed to

describe the differences between an object with an approximate dimension of 2.37

and another with a dimension of 2.45. The pictures are much more descriptive.

One use of approximate fractal dimension is to describe an object's

relative roughness. It is beneficial to view a fractal dimension as degrees of

roughness between the standard three dimensions. If the dimension is between 1

and 2 then the object should be a very irregular curve. If the dimension of that

curve approaches 1 then the curve is probably not very rough and would lack any

interesting diversion from an ordinary plane curve. If the dimension of the curve

approaches 2 then the curve becomes like a plane or filled polygon and again

lacks appeal. The most interesting fractal curves are those which demonstrate

dimensions nearer the center of the scale between the standard Euclidean

dimensions. A similar argument can be made for solid objects with dimensions

between 2 and 3.

This is not to say that fractal geometry is not a powerful tool for the

graphics programmer. The evidence of the power of fractals to model objects of

considerable complexity is clearly demonstrated. To date, this power has not

been matched by other standard methods.

The graphics programmer should not concern himself greatly with the

dimension that is demonstrated at different levels of object construction. He must

concern himself with the techniques of construction and the realism that is

achieved.

7
Mandelbrot's method for estimating the Hausdorff- Be.sicovitch dimension for non-random
sets built through self similar shapes will be introduced in section C after the hoch curve i>
described.

23
C. FRACTAL CURVES AND SETS
1. Definition of Fractal Sets

We take our definition of fractal sets from [Ref l:Chap 3 and pp. 361].

Definition:

A fractal set is a set for which the Hausdorff Besicovitch dimension strictly
exceeds the topological dimension.

As we have established and is emphasized by Mandelbrot, this definition

is not very useful. The definitions of topological and Hausdorff dimensions are

very involved. It is a gargantuan effort to prove that a set has a topological or

Hausdorff dimension (if one desires complete rigor, typically the topological

dimension is derived by the least parameter approach (section 2)). When using

fractal functions then, it is practical to assume that because the functional

method you use is fractal-like that the dimension is fractal.

The functional techniques to be introduced have a certain methodology

that creates fractal sets with a behavior that is disciplined and predictable. The
assumption is, since these methodologies are well behaved, that any set created

by these methods (with some careful restrictions) will itself be a fractal set.

WeT
are left with a practical methodology whereby we discover fractal

functional methods, prove that the set created is fractal, characterize the fractal

part of the functional method (carefully) and then enshrine that method as a

8
fractal method. If one's purpose is a practical application of fractal techniques

and not a rigorous mathematical investigation then this is a reasonable and

practical approach. This approach is taken in the remainder of this study.

2. Constructing Fractal Sets

In order to describe the construction of the Koch curve, it is necessary to

present terminology introduced by Mandelbrot [Ref. l:pp. 34-35].

We use a geometrical shape (at first a straight line) and call this shape

an INITIATOR. We create another shape that is constructed with shapes

similar to the initiator and call this a GENERATOR. We define a sequence of

8
A working model or equation where you are only concerned about the behavior and not the
exact mathematical properties.

24
transformations upon all current initiators (suppose there are m such initiators)

in the construction by applying the generator to all initiators. This creates a new

construction that consists of m x r (where r is the number of distinct parts of the

generator) sides where each side is a shape that is similar to the initiator. The
next step is to again apply the generator to all initiators.

This recursive definition has no terminating event but is continued ad

infinitum. This functional process is well suited for recursion because the

application of the generator to the initiator is constant with respect to method

and varies only to scale. It is also well suited for parallel processing (in the

computer science sense) because each application of a generator on all current

initiators is independent.

These concepts are probably confusing at this point and were especially

difficult to visualize when they were first envisioned because the authors had few

tools beyond mental imagery to convey their point. This is probably why they

were largely ignored for 70 years. It is much easier to visualize these functions

when they are shown on a computer graphics display.

3. The Koch curve.

The mathematicians Koch and Cantor developed functions which

attempted to challenge the mathematical models of continuity and

differentiability. These equations were developed during the great debate on set

theory and were used by Cantor in arguing for his theory. These functions were

like none before, using constructions which played upon natural geometric

constructs but when combined with the power of infinite recursion became sets

which defied intuition. It was not until much later that mathematicians were able

to reconcile these functions with algebra and set theory. The function to be

introduced has been proven to have a Hausdorff-Besicoritch dimension which

exceeds its topological dimension .

The Koch curve is a very beautiful curve that at first gives the observer

the impression of a snowflake or a coastline (Fig 2.1). Mandelbrot uses this type

of construction (with variation — i.e. randomness or less behaved generators) to

draw coastlines that look very realistic. To understand this construction is to

Thus fractals do exist and it is possible to rigorously prove it to be so.

25
understand one method of obtaining a fractal set from a well defined non-fractal

set (R 2 ) using non-random techniques. This method of construction (in the

general sense) is very powerful and is used throughout this thesis.

To construct the Koch curve, we use three initiators (line segments) of

10
equal length and join them to form an equilateral triangle . To construct the

generator we use four line segments that are each — the length of the initiators

and apply these to each initiator (Fig 2.1). This yields a new geometric figure

with 12 sides versus the original 3 and a total perimeter length of 4 units of

length versus the original 3 units.

Figure 2.1 demonstrates the first and second recursive iterations of

building the Koch curve. Observe how the progression develops to yield the final

figure in Figure 2.1. Imagine this progression occurring indefinitely.

- With each iteration of applying the generators the total perimeter length

increases by — over the previous perimeter.

- The length of the curve begins to increase without bound even though the
length of the initiator decreases to an infinitely small length. Hence the
curve's length unbounded with no point
is intersecting but yet is contained
in a small bounded two dimensional area.
- The points of the curve are by construction only the end-points of each
initiator and each point is clearly distinct from the other ( no two points are
connected).
- Although each point is distinct at any one level of the curve construction, it

can be proven that the curve when viewed in the limit is continuous at every
point.

- That due to the above qualities the curve is not differentiable at ANY point.
It is important to realize that the endpoints of the lines (initiators) are the only

points of the curve. The line only serves as a vehicle by which the points may be

easily determined. The exact same set could be built using 180° arcs as initiators.

An algorithm and computer graphics program for the construction of the

Koch curve is presented in Chapter 4 and Appendix A.

The choice of an equilateral triangle was arbitrary. We could have chosen any shape as
long as it was made up of Initiators and avoided intersecting lines during recursion.

26
 :

INITIATOR

GENERATOR N = 4

r = -

THE KOCH CURVE

Recursive term- Output medium:

ination distance Laser printer
is .05 inch. Resolution
300 dots/inch
DT = 1 , HB = 1.26

Figure 2.1 The Koch Curve

27
 4. Mandelbrot's Dimension Approximation Function

In the above section, one functional construction technique for building a

fractal set has been introduced. It is possible to approximate the dimension of a

fractal curve that is built using these constructs. In [Ref. l:pp 56-57], Mandelbrot

introduces a function that is based on the similarity properties of the above

technique. This function has a real exponent D that is the approximate

11
dimension of the fractal set .

Consider a method of paving (covering?) a Euclidean shape. Divide a line

segment into N segments with each segment a part of the original segment such

that the sum of the lengths of the N segments equals the length of the original

segment. It follows then, that the sum of the ratios of the divided segments
/ . . \
Divided segment
lengths to the original segment length e.g. rs = —— ;
s
must
Original segment

equal 1.

N
£r s
= 1

s=l

We know that the dimension of a line is equal to 1. If we raise each of the above

ratios to the power D (where D = 1) the equivalence still holds.

Eh =1

Let's allow the Koch function to assume a similar dimensional relationship but
treat D as a real valued unknown. Refer once again to Figure 2.1. Notice that

the length of the four line segments that make up the generator have a length

ratio of — to the initiator. Call this ratio r m . Notice that the generator is made

up of four initiator shapes. Call this number N.

s M = 1

11
CAUTION: This technique does not necessarily apply to other fractal functional
methods.

28
 )

Substituting and Solving for D:

D
|rJ -

/ \ D
1
4 x = 1
3

D = -j^- % 1.2618
log3

Which is equal to the Hausdorff-Besicovitch dimension of the Koch curve 12 .

This is not a proof of a general equation for the fractal (Hausdorff)

dimension of a self similar fractal set but implies that a general dimension

generating function is possible:

G(D) = £ (r m
m=l * j
'

Where N = number of sides of the generator.

Where Rm = ratio of side m to the initiator.

Mandelbrot claims that experimental evidence suggests that this equation holds

whenever this functional method is used.

When each segment of the generator is a fixed ratio to the initiator (as is

the case with the Koch curve) then the solution to this equation is trivial:

lo 8(N)
D
l°g(—

Functional Characterization of the Koch Method

A complete and rigorous inductive definition of most fractal functions

can stretch the notational capabilities of the symbolic aspects of the inductive

and functional methods. It is thus generally impractical to use these methods

I2
The Koch curve was proven to have a Hausdorff-Besicovitch dimension equal
O £J 4
I

~ 1.2618 and a topological dimension equal to 1, Hausdorff Dimension und au>seres \l iss
log3

29
except in a verbose and non-rigorous manner. It can be insightful however, to

dissect the beast (once!) and hopefully gain further intuitive insight.

To simplify the process, we define a Koch half-line as in the above fractal

but with a single line segment as the Initiator (versus the equilateral triangle).

We consider the line interval of [(0,0), (1,0)] within the set R 2

. This restricts the

fractal shape that is drawn to a partial function on: [0,1] X [0,1] -* [0,1] X [0,1].

Using the inductive process to define the essence of the fractal

sequencing, we have the following definition:

(Basis)

Step

n = {(o.o,o.o), (i.o,o.o)}

(Induction)

Step k

Label the ordered set of 2-tuple points from" Step k-1 as:

n*-i = {P k -\P k -\P k -\,..,P k -\ }

where
n = 4(*-l) + i

Determine the new set of ordered points for Step k as:

n* = {P* P* 2 P* 3 ,...,P* m
1 }

where
m = 4* +1
where
P*! _
iJfc r>k-l
= P K ~\

Pk 5 = P k '\
Pk o = P k ~\

pk _ pk - 1

30
 where

(p\,p* 7
,p* 8)
= P k ~\
[
2
R ]p*
3
'

~ k-\
C* 10/* 11/* 12) - Pk X
:
3R4

I pk pk pk ok -1 *-]
V-* m-3, r m-2, r m - l)=
1 ^ n-1 n-l^nj

and where the relation ,-iR. is the geometrical relationship between the two

end points of the initiator line segment from step k-1 and the five points of the

generator that compose the ordered subset for step k as above. For purposes of
brevity, the full functional definition for this geometrical relationship is explained

in detail in Chapter 4.

(Extremal)

No point is an element of Q unless it can be shown to be so in a finite

number of applications of clauses 1 and 2 above.

It is thus possible (but tedious) to rigorously characterize the Koch
fractal set within the well defined constructs of induction and the functional

technique.

6. Another Fractal Set, Cantor's Dust

Cantor developed a function 13 which used the same functional technique

as the Koch curve but with a reverse twist. Cantor's function takes an initiator

(the unit interval [0,1]) and dissolves it into a discontinuous set which is as rich

in points as the interval [0,1] but contains no interval itself [Ref. 6:pp. 22-23].

The points of the set are all distinct but the set has the same cardinality as the

unit interval. The best description for this set is that it resembles a dust.

is
Also referred to as Cantor's discontinuum or Cantor's triadic set.

31
 Cantor's Dust is difficult to demonstrate on a graphics display because it

quickly dissolves below the resolution of the display. Refer to Figure 2.2 to

visualize the initiator and the generator. The initiator is a line interval [a ,6 ]

where ^ a < b ^1 and the generator is two intervals each — the length of

the initiator such that interval 1 is [a , a + (b -a) x — and interval 2 is

[6- (b
V
-a) x — ,6]. After the initial application of the generator to the unit
'
3

1 2
interval there will be two intervals in the current construction [0.
— ] and [
— ,1].

The second iteration of the recursive routine will yield four intervals

2 1 I t 8
[0.-1, [-,-].
1 J l J
[-,-] andJ
f-,ll.
l J

9 9 3 3 9 9

Every initiator that is created by the generator has an infinite sequence

that is begun at the next application of the generator. This causes a series of

convergent sequences to each end point of each initiator. Thus each initiator

= 2
CANTOR'S DUST
INITIATOR £
GENERATOR * = °
5 1
i

• • • Initiating Structure: D
T =
The Interval [0,1]
HB = -6309
1st Iteration Initiator Length ^

2nd Iteration . Initiator Length -

3rd Iteration Initiator Length ^_

4th Iteration Initiator Length ,5-

o1

Figure 2.2 Cantor's Dust

32
spawns a convergent sequence toward its endpoints. But for each initiator there

are also two spawned sequences toward the center. It is possible to imagine this

as an infinitely dividing organism which leaves behind four points (eggs) each

time it divides and then each egg itself replicates • •

ad infinitum.

This functional method of dissolxnng a line is very powerful as a tool in

computer graphics because it can be used to cause many special effects from

ordinary objects. Mandelbrot has used variations of this method to create images

of star systems for example.

The topological dimension of Cantor's dust is and by Mandelbrot's

dimension generating function we have the fractal (Hausdorff) dimension equal

to:

N / m
G(D) = £ rj = 1

m=l V
'

Where N = 2

Where Rm = -
3

Substituting and solving

D
'i
2 x = 1
3
v /

log 2
D = -^£_=. « .6309
log 3

Cantor's dust was proven to have a Hausdorff-Besicovitch dimension of

-i^Jl ~ 6 3 9 [Ref. l:pp. 77-78].

log 3

7. A Note on the Concept of Similarity and Fractals Sets

The use of the mathematical concept of similarity continually shows up

in the investigation of fractal functions. This can be expected because of the type

of functional building tools that are used for most of these sets. The relationship

of the generator to the initiator has similarity built in. Thus Mandelbrot is able

to use similarity properties in formulating a dimension generating function and in

making claims about the set's inherent structure.

33
 The functional method invented by Koch and Cantor is but one method

of determining a fractal set. Similarity fractals may eventually be grouped into a

class of fractals (important but restricted). This type of functional method

provides a vehicle for the creation of disciplined fractal sets but makes you

wonder about the rest of the space that fills the gaps between the standard

Euclidean shapes 14 the self-similar fractals and the infinitude of

, R 1
.

Much like the transcendental numbers (n and e for example ) in all their uncountably rich
expanse fill the gaps between the familiar (and well behaved) sets in R.

34
 III. THE IMPLEMENTATION OF FRACTALS IN COMPUTER GRAPHICS

This chapter introduces the reader to the practical aspects of implementing

Fractals on a two dimensional computer graphics display. We are forced to

confront the issues of economy of scale between the infinite fractal function and

the finite computational environment of the computer. L'nderstanding these

compromises is a necessary bridge to successful fractal programming.

A. THE IMPLEMENTATION PROBLEM

1. Infinite Recursion, Stacks and Data sets

It is naive to view the computer as a truly infinite abstract machine

which is capable of any binary computation of any length. Infinity is a concept

that when applied to physical objects quickly breaks down as soon as that object

is bounded in any way.

It is possible to model infinite behavior in computers though mathematics

(automata theory) and gain useful insight into possible capabilities of the

computer (the use of push down automata in compilers for instance). But in

order for automata theory to make assertions it is often necessary to make the

assumption that the automata (computer or an abstract machine) is in fact

infinite. When the assumption of infinity is made, there are many powerful

mathematical tools which can be brought to bear upon non-intuitive abstract

13
problems that would otherwise be functionally intractable if the automata was

considered to have an arbitrarily large but bounded space. The question arises,

can we consider the implementation of these constructs (insights) as valid? The

answer is yes. but only in the context of some bounded space (for instance, the

maximum stack space for a push down compiler).

The question of how many valid programs can be recognized by such a

compiler is usually too difficult to determine. The compiler's solution to the

problem is to use a passive sensor to detect stack overflow and notify the user

15
Functionally in the sense that the problem may be decidable but the solution space is so
large and undefined that it could not be determined in a reasonable amount of time.

35
that his problem is to large for the current stack space. Although the compiler is

theoretically a recognizer for an infinite set of programs, the finiteness of the

computer is the grounding factor.

A similar paradigm exists for implementing infinite functions like fractals

in computers. Fractal functions use recursion, randomness, massive floating point

computations and large amounts of primary memory. In order to use the fractal

function productively, we must manage the methods we use and produce a finite

approximation of the fractal set we create. This can be done by a passive or

active means.

2. The Bounded Stack Limits Recursion

The Koch-like fractal functional method is by definition a recurrence

equation. All fractal methods introduced in this thesis have a similar recurrence

definition. Thus it is impossible to avoid the use of recursion in the production of

16
fractal graphics .

Recursion on most Von Neuman type computers is implemented on the

system stack (which may be hardware (fixed) or software (in primary memory
and therefore expandable)). Such a stack is always bounded by a fixed upper

limit of allocatable memory space. The formal definition of the Koch-like fractal

method has no recursive terminating event. In order to use these methods, we

must determine the precision that we require and develop a termination event to

signal the beginning of the recursive ascent when this precision is met.

Failure to manage the recursive descent inevitably causes the exhaustion

17
of the system stack which stores the program state on each recursive call . The
programmer should be careful to keep his local variables at a minimum in the

recursive subroutine. By doing so, the total data stored during each recursive

call is reduced.

16
Why should you try ? The recursive part of fractals functions is the mathematically beau-
tiful aspect which makes them so eloquent and powerful.
This could be the recursive termination event if you are not careful. It might be a good
idea to try this experimentally on your computer to determine this maximum recursive descent
distance.

36
 3. The Computer Graph i cs Set Paradigm
The paradigm that is used for displaying objects on the raster graphics

screen is inherently discrete and two-dimensional. A typical system consists of

primitives that allow the user a view of his modeling world as a largely

unrestricted three dimensional space. These primitives are limited in their power

and are usually based on the Euclidean geometry (lines, points, polygons etc.).

The user is required to supply viewing parameters that define a limited three

dimensional viewing space. These commands are processed by the graphics

system which projects the objects contained in the viewing space onto a two-

dimensional space that is the display screen. The display screen can be divided

into discrete entities called pixels, each of which represents one point on the

screen. The number of pixels defines the resolution of the screen.

Given the graphics paradigm, it is natural to model objects that are

defined on the display as a collection of discrete points. When viewed on the

screen, these pictures can accurately model objects of considerable complexity

[Ref. l]. This view is a departure from the normal view of the graphics

application programmer. The application programmer views the objects he

creates through a collection of Euclidean geometric primitives (lines polygons

18
etc.) that are abstracted from the actual display . Both views are useful models

in describing objects but for the fractal graphics programmer the former is the

more powerful because it neatly maps to the fractal method of object

construction.

For this study, we view the world and the objects it contains as

collections (sets) of discrete three dimensional points. When an object is built

through fractal techniques, the length of the pixel is a natural termination point

for the fractal recursive process. This yields an attractive bijective mapping from

our object to the display screen.

18
All computer graphics is an abstraction; realism is achieved through deception of the eye
by a very small collection of colored points. Our goal then, is to find the most efficient method of
defining those points.

37
 4. Summary
From our previous discussion, it should be obvious that fractal sets do

not exist in our computer. The sets that we create are finite approximations of

the actual set. The same is true of the fractal pictures that we create and to

which we ascribe fractal dimensions. It is an illusion of the eye that we create by

taking the fractal computations below the resolution of the screen.

B. MAPPING FRACTALS TO A BOUNDED SPACE

In view of the graphics set paradigm, it behooves us to develop a completely

bounded set space which has a one to one correspondence between the points

computed by the fractal equations and the entities of the screen. This abstraction

is appropriate because these entities called pixels (which are nothing more than

colored points) are the only components of the picture 19 . The only reason for not

using this methodology before is that it was functionally difficult to compute

(without some simplifying abstractions) the large number of pixels in the display

set.

1. Fractal Recursive Termination Event

Most graphics software packages provide a means by which the user can
define the space on the display screen onto which he wishes to map. This is done

by providing to the graphics system, viewing parameters, which define a bounded

three dimensional space which is then projected onto the two dimensional screen.

Many mappings are possible within the above constructs but for the

fractal programmer it is convenient to establish a mapping that provides for a

one to one (or bijective) relationship between the 3 pixel set and the real-valued
coordinate space in which we compute a fractal object. We develop this bijective

relationship by defining a fractal part to be determined when the distance

between generating points (in real coordinates) is less than the distance of one

pixel in the mapping of real to screen coordinates (SCR), Figure 3.1. This

mapping is explicit when the window, viewport and z- clipping are defined. This

mapping limits the number of pixels that can be associated with a given fractal

initiator to a one to one mapping of fractal parts to the pixel set. The size of the

In fact, this is the way the Euclidean graphics abstracts such as a line are mapped to the
screen. At some point in the display process this mapping must take place.

38
pixel is a cubic space which clumps together the infinite fractal space into a

discrete number of cubic spaces equal to the number of pixels in the pixel set.

The fractal programmer should choose his viewport and window carefully

so as to elicit the most attractive mapping possible between his object and the
screen. The viewport (SCR) rectangle should be defined such that the ratio of

the X distance (horizontal) to the Y distance (vertical) is equal to the ratio of the

X-Y distances of his (real-valued) window. The Z-distance (front and back
clipping plane distance) should be equal to the distance in the window
coordinates that defines the length of the pixel multiplied by the desired Z-depth.

Figure 3.2. Formally:

X SCR WLD
SCR Y WLD

Recursive Termination Event

Recursive Iteration i Recursive Iteration i At
The POINT PLOT Routine

z
2
The Initiator Length Is Less The Pixels Were Mapped To
Than The PIXELSIZE, So Map Their Corresponding Position
The Points To The Pixel Set. In The Pixel Set.

Figure 3.1. The Fractal Recursive Termination Event

39
 1 ; —

Graphics Display
Terminal The Pixel Set
'User Defined Viewport

'Graphics Display Space

1 1 i
'
I : ,
'
'
'
1

1 I
•
i
'
'

: ! i ' : l
'
*-v~-
h :
'

;
:,
M
i
;
i
i i

User Defined Window \ I

i

i
1
i
'

i
'
.
i i
:
i < i i
'
!
'
'

:
i ! I
i
!

«e
|

In World Coord. System \

-Ur-
'

'J.
Y Distance
II'' w *~
—
'ii
in WLD

=
' '

^ii i
'iifr?
T iT
-» t
H
I I I :
i
i i i i
:

..-K ;. -
.

5^E
3±Z t '
'

I I I I
i
' 1 I

i ' i
| I I I

X Distance
WLD

The Cubic Pixel Space

The Pixel Size Is Defined By The Z Depth In Pixels Is Defined By

XwLD Z-Depth^D
PIXEL SIZE = Z-Depth pij
SCR PIXEL SIZE

Figure 3.2. Normalizing the Pixel Space

40
If the above relationship holds then the size of the pixel in world coordinates
equal to the X (or Y) world distance divided by the number of pixels, Figure 3.2.

PIX SIZE =
XsCR

The front and back Z-clipping planes should be established such that the distance
from front to back is equal to the desired Z depth (expressed in pixels) multiplied

by the pixel size (in world coordinates), Figure 3.2.

If the above viewing relationships hold, then the inclusion of a simple

check:

IF (InitiatorJDistance < PIX_S1ZE) THEN

return;
ELSE
continue;

terminates the current recursive descent and begins backtracking.

2. Memory Requirements
The amount of main memory required by the fractal programmer is

directly proportional to the size of the pixel space that is being used to display a

fractal figure and the sophistication of the desired display,

a. Data Locality During Recursive Descent

If the fractal in question is a 1-2 dimensional display that is

displayed as a collection of points on a two dimensional plane then the

requirement to store the 3-tuple points is eliminated ((x .y ,z ) coordinates which

represent a point in R 3
). The points can be computed by the function and
displayed on the screen in the immediate mode and then destroyed. The only

requirement for memory is the data that is germane to the program ("globals"
and thus fixed at run time) and the amount of data pushed onto the run-time

stack for the recursive calls up to the point of the deepest recursive level. This

requirement is minimal and does not present any significant limitations.

An example of such a fractal is the Koch curve. The algorithm

presented in chapter 4 uses the data that is local to the subroutine that is

computing the current generator with the input coordinates of the initiator. The

41
inductive nature of the fractal method provides data independence so that the

programmer does not have to have available to his subroutine all points

computed thus far.

b. Memory Requirements for the Fractal Set Paradigm

If the programmer requires a more sophisticated display that includes
a requirement for hidden surface removal and (or) lighting enhancement then the

fractal method becomes just a part of the overall display process. The
programmer has the fractal process compute the points of the picture and these
points are stored in a memory structure. This structure is processed by a hidden

surface algorithm and a lighting algorithm and is then projected onto the screen.

The ideal Fractal computer would have an exhaustive memory and

unlimited computing power to be able to allow us to store the entire pixel set in

memory. This however, can stretch the capabilities of most present day

computers. For example; if the pixel set is defined such that its dimensions are

1000x1000x1000 (which is certainly reasonable), the amount of storage required

20
would be an array of 10 9 x 24 bits or 24 billion bits.

Current techniques exist where this storage can be minimized. A Z-

buffer array can be used to store the current Z coordinate of the forward most

displayed pixel. By using this method the fractal computation can be made and

then a hidden surface computation can be immediately invoked to check the

visibility of the point by checking it against any other point in the same position

on the X-Y plane. This technique reduces the above space requirements to

6 6
10 x 24 bits + (Z- buffer = 10 x 32 bits) 21 or 56 million bits.

A fractal programmer must manage his memory resources carefully

in order to maximize the computational resources available to him. This usually

requires a tradeoff between efficiency, realism and the memory space utilized.

The algorithms for hidden surface and lighting are well covered in the literature

and although they are integral aspects of the fractal realism issue they are not

germane to this study. Throughout the remainder of this study, we are not

" 24 bits for a machine assuming the RGB color system with 24 bit planes. Each color; red,
green and blue would then have 8 bits of precision.
1
Assuming the Z coordinate is stored as a 32 bit floating point number: this number can be
further reduced if we restrict the Z precision.

42
 l

concern is of cr. ng the --ume *.

have the full three dimensional •

3. Concurrency

A question arises is it possible to |

- ..-. Real Tim". We
have enough knowledge at this juncl in I liscuss the "or concurrent

operations during the Brad iescen: 7 gies

that can used if concurrent operations are *ompu:er We could

assign a process to each of the initial .*•

r ts d -.red

precision. We could also h. . : rocessor U "•e am- .ich

at any level of the i irshn mil pri-ce^- fcc

_ erator to the local initiator"

The local data inde: taJ rec orsr a allows such

i sj awning because the application :'':.• gena I - "o an in:" "

i - :<al :c

initiator. This applicati u is i res ill : the :.'.'. -

dT the K-l -

thai led up to the K :r sterf ~'r.~ mathematical process ::" ap; .y.r.g :':.< - -

the K level needs only the initiator data from the K— - step and reqoin

knowledge of the entire :"rac*a. sel Tne ••- 1 indue* ion - "hen

computationa..; .ized. The number I as bel n :he en*

determined is not _: : than the number of computations required to app

generator times the ~

"
: . .:.. recursive - . \:. : :.
-
:or

recursive termination event. A short algorithm describing the processing . **s

of the second meth 1 is shown in Figure 3

" To beat the hur -rid ach.-: graphics illusion, typkalh by completing

all computations and screen I O prior to the — •• -

"or a SOHZ display) of most

raster graphics i

- of course would require an enormous n rower thai does no* corr-

M in reality the spawi r - -Id have to be bounded
- in some -
 main()

LOAD INITIATOR AND GENERATOR DATA

concurrent;
generate (INITIATOR l)

generate (INITIATOR 2)
generate (INITIATOR 3)

generate(INITIATOR n);
end concurrent;
end main;

generate(INITIATOR I)

BUILD GENERATOR FROM INITIATOR

if (Distance < PIX_size) then
PLOT POINT
return;
endif

concurrent;
generate (INITIATOR 1)
generate(INITIATOR 2)
generate(INITIATOR 3)

generate(INITIATOR m);
end concurrent;
end generate;

Figure 3.3. An Algorithm Which Utilizes Concurrent Operations.

44
 IV. FRACTAL COMPUTATION IN R 2

The algorithm introduced in the following section is capable of computing a

very broad class of Koch-like fractal curves within R 2

. It provides the fractal

graphics programmer with the basic tools for fractal computation and an
algorithmic template that can be used for many applications. The graphics
programmer needs to fully understand fractal programming within R" before he

attempts the more complicated concepts of fractal terrain modeling in R3 .

A. THE GEOMETRY OF INITIATOR -» GENERATOR

The geometric relationship between the Koch curve initiator and generator

can be described through a very simple set of data that captures the essence of

the Koch curve. The method introduced also allows us to vary the data which

defines the generator and compute many different fractal shapes .

The general strategy for computing generator points from a set of initiator

points is to determine two lines that intersect at an unknown generator point.

Figure 4.1. a. The unknown generator point is defined by constant relationships

between the initiator and generator. The first line is the perpendicular from the

generator point to the initiator line. The second line is formed between the first

point of the initiator and the unknown generator point. All data to compute the
line equations are derivable from the two endpoints of the initiator or from

constant initiator/generator ratios. For simplicity's sake, we ignore the divide- by-

zero problems that are encountered when any of the lines are parallel to the X or

Y axis. These situations only simplify the computations and their solution is

demonstrated in Appendix A.
If we label the endpoints of the initiator as P =(X1 1
,Yj) and P : =(X : .Y : )

(Figure 4.1. a) then the slope of the initiator line is:

Y, - Y,
Slope, init =
X2 - X,

We must be careful in our terminology because the relationship described can draw shapes
that do not avoid self intersection and thus must be considered quasi-fract al.

45
 The slope of the perpendicular intercept line is the negative inverse of the

slope of the initiator, hence the slope of this line is:

Slope, perp =
Slope. init

The intercept point on the initiator can be determined by using a constant

distance ratio as shown by the following equations:

+ X2
X.perp = —Xj (Generator. ratio. constant x )

1 + Generator. ratio, constant

Fj + [Generator. ratio. constant x Y2 )

Y.perp =
1 + Generator. ratio. constant

The value of " Generator. ratio. constant" is a fixed constant (although it might be

interesting to randomize it) where you determine at what point the generator
intercept point intersects the initiator and express it as demonstrated in Figure

4.1.b. This ratio can be determined graphically, through hand calculation or via

an interactive automated means.

With the slope and a point of the line (X.perp, Y.perp) we can determine the

line equation for the perpendicular line by determining the Y intercept:

Y. intercept. perp = Y.perp - (Slope. perp x X.perp)

This yields the line equation for the perpendicular line:

Y = (Slope. perp x X) + Y. intercept. perp

To determine the generator line equation for the line segment which connects

the first point of the initiator and the unknown generator point, we need
constant information about the angle between the initiator and this line. This

angle is always constant with respect to the initiator and like the ratio

information above, it is recorded as a constant at run time (see Figure 4.1.c),

note that the angle may be positive or negative.

46
 X Yg.j>Jif*+ Unknown
( gmn '
Generator Point
Generator Lin<
(x i >
x*>*
Initiator

(X >
Y J :

•Perpendicular Intercept

Fig 4.1. a Intersecting Lines Determine a Generator Point.

'unk

Initiator Generator Intercept Point

Gi
dist. a. dist . b
H
PiGi a.
Ratio-constant =
Gi P2 b

Fig 4.1.b. The Generator Ratio Constant

Point with positive angle Point with negative angle

.' G,

/ e

\ e

Fig 4. I.e. The Angle '\/*a

Figure 4.1 Building the Generator with Intersecting Lines

47
 We record the data about the angle as the tangent of the angle. With this

information, the slope of the generator line can be determined with the following

equation:

01 Tan 6 + Slope. init

Slope. gen =
(1 — Tan 6 ) x Slope. init

The Y intercept for the generator line can now be determined:

Y. intercept. gen = Y — l
(Slope. gen x X x)

This yields the line equation for the generator line:

Y = (Slope. gen x X) + Y. intercept. gen

The Cartesian points of the unknown generator point can now be determined
by intersecting the two line equations and solving for X (Figure 4.1. a) then

substituting X into one of the line equations and solving for Y . The
equations follow:

v
A8 =
,. .

Y. intercept. perp —
Y. intercept. gen
Slope. gen — Slope. perp

Y gen = (Slope. gen x X gen ) + Y. intercept. gen

The constant data for the Koch curve that corresponds to this geometric

method is illustrated in Figure 4.2.

There are many different ways to build a generator given an initiator using

standard geometric constructs but the method introduced allows experimentation

with the Koch function to discover new shapes. By varying the tangent of 6 and

the fixed ratio of Figure 4.1 we can describe new generator constructions. These

new constructions can be used in the same algorithms that compute the Koch
curve. We initiate the recursion on a generating structure built of line segment
25
initiators and allow the recursion to progress until a terminating event , creating

many diverse shapes. Figures 4.6 through 4.9 demonstrate some of these shapes.

The length of a pixel or an arbitrary line length.

48
 The Koch

Point
_W_
Point
Generator

i
1 2 Point 3
i
i

i i

i i

P
>4-J P Pl
i

G P
2 3 8 3

Ratio = O.S Ratio l.O Ratio s 2.0

0= o.o 0= 0.4718 radi. 0- O.O

Figure 4.2 Ratio Constants for the Koch Curve Generator.

B. THE MID-POINT DISPLACEMENT TECHNIQUE

Mandelbrot [Ref. l:pp. 43 and pp. 233-234] uses an alternate technique to

draw the Koch-like curves that is equally valid. He calls his technique mid-point

displacement because he determines the generator via fixed relationships that

displace a point from the mid-point of the initiator. This method uses many of

the same geometric relationships that are used above but provides a different

progression to the method of building the generator. This new view allows us to

look at the relationship between the initiator and generator in a slightly different

light. By so doing we are provided new insight as to how we might alter the

relationship to create new images and sets.

The mid-point displacement technique can be useful for two other reasons. It

is the best known method for fractal set building and because of this, it facilitates

communication between fractal programmers. Most of the terrain models that

have been developed use the mid-point technique. It also provides an easier

method to avoid line intersection.

49
 The Midpoint Displacement
Technique
,' i

Figure 4. 3. a First Midpoint Displacement

Figure 4.3.b Second Midpoint Displacement

THE KOCH GENERATOR

V--- ^\
Figure 4.3.c Third Midpoint Displacement.

Figure 4.3. The Koch Generator Using Midpoint Displacement.

The mid-point displacement method is demonstrated in Figure 4.3. The

method progresses by taking the initial initiator and applying the first midpoint

displacement. This yields the figure demonstrated in Figure 4. 3. a. The next step

performs a mid- point displacement on the left initiator created by step 1. This

yields the figure demonstrated in Figure 4.3.b. The third and final step is to

replace the right initiator created by step 1 as demonstrated in Figure 4.3.c.

50
 a

The inversion of the direction of the mid-point application (in Figure 4. 3.

the displacement is above the initial initiator where in Figure 4.3.b it is below)

can be accomplished with a single computational procedure. We only need to

invert the position of point 1 and point 2 in the parameters of that procedure.

The parameter inversion changes the orientation of the initiator in space with

respect to the computation of the midpoint displacement. The procedure blindly

computes the mid-point displacement relative to a fixed relationship to the

initiator input points. This operation implicitly defines the orientation of the

generator in space.

The geometry for computing the midpoint given two initiator points can be

computed in precisely the same manner as the intersecting line algorithm (The

Koch midpoint ratios are 1.0 and the angle is .437 radians). An alternative

method is to use the equations of a line (or a plane) normal. The second method
(utilizing the normal) provides a geometric relationship which is intuitively

appealing. Its appeal comes from the desire to modify the length of the

displacement relative to the initiator with a random scaling factor.

The mid-point displacement technique has some advantages over the line

intersection algorithm. Random modification of the length of the displacement

along the normal (from the computed generator point to the initiator using the

line intersection algorithm) is not intuitively appealing. It requires a translation

of the desired displacement into an angle (the angle 9 (Figure 4.1.c) between the

initiator line and the unknown generator intercept line). The control of that

angle is less intuitive than the control of a displacement length. The geometry
for mid-point displacement using the initiator normal is introduced in Chapter 5.

C. A KOCH-LIKE FRACTAL ALGORITHM

Implementing the above function is a relatively easy process that is

demonstrated in this section and Appendix A via a gradual unwrapping of the

layers of complexity that are required to successfully implement the algorithm.

The algorithm roughly follows the template used in Chapter 3 to demonstrate

concurrent processing. A C-like language is used for the algorithms.

The first algorithm (Figure 4.4) is a template that delineates the basic

processing steps. This recursive process is typical of fractal functions and can be

51
used as a template for many fractal programs. The second algorithm (Figure 4.5)

is an expansion of the first and demonstrates the replacement of a given initiator

using the line intersection method.

Appendix A is a complete Fractal program. This program was used to

produce the data for Figure 2.1 and Figures 4.6 through 4.9. This program

demonstrates the precautions that must be taken to avoid divide- by- zero when

lines are parallel to the X or Y axis.

main()
{

Load Initiator Coordinates;

Load Generator Relationship Values;

for 1=1; I<= Number of Initiators; 1=1+1;

{
generate(X(I) 1 ,Y(I) 1 ,X(I) 2 ,Y(I) 2 )
}

generate(X 1 ,Y 1 ,X 2 ,Y 2 )
{

Determine Distance Between the Endpoints of the Initiator

if (DIST < Pixel.length) return;

Replace Initiator with the Computed Generator;

Load the New Initiator Data into Local Generator Array

for J=l; J<= Number of Generator Segments; J=J+1;

{
generate(Xgen(J)!. Ygenp^, Xgen(J) 2 Ygen(J) 2 );
,

Figure 4.4. High-level View of the Koch Algorithm.

52
main()
(

Load Initiator Coordinates;

Load Generator Relationship Values;

for 1=1; I<= Number of Initiators; 1=1+1;

generate(X ,Y ,X 2 ,Y 2
1 1 )

/* Determine Distance Between the Endpoints of the Initiator '

/
DIST = sqrt((X 2 -X 1
)**2 + (Y 2 1
)**2); -Y
/* // Distance is less than Pixel Length; Plot and Return */
if (DIST < Pixel. length)

{
plot.point(); /* Your Graphics Point Plotting Routine */

return; /* Point 1 and 2 plot the same pixel */

>

I* Load The Endpoints of the Initiator into the Generator Array '/

Generator. X[0] = Xp
Generator. Y[0] = Y,;
Generator. X[Number.of.generator. points + l] = X: ;

Generator. Y[Number.of.generator. points + l] = Y2 ;

/* Determine Slope of the Initiator */

Slope. init = (Y 2 -Y,) / (Xj-XJ;

- continued-

Figure 4.5. Detailed \'iew of the Koch Algorithm.

53
/* Calculate Unknown Generator Point via Intersecting Lines
the */
for J=l; J<= Number.of.generator. points; J=J+1;
{

/* Determine the Generator Point Intercept on the Initiator */

X.perp = (Xj + (G enerator .ratio. constant [J] * X 2 )) /
(1+ Generator. ratio. constant [J]);
Y.perp = (Y l + [Generator. ratio. constant [J] * Y 2 )) /
(1+ Generator. ratio. constant [J]);

/* Determine the Slope of the Perpendicular Line */

Slope. perp = (-1 / Slope. init);

/* Determine the Y -intercept of the Perpendicular Line */

Y. intercept. perp = Y.perp - (Slope. perp * X.perp);

/* Determine the Slope of the Generator Line */

Slope. gen = Generator. tan. theda [J] + Slope. init) /
(1- Generator. tan. theda [J] * Slope. init);

/* Determine the Y - intercept of the Generator Line */

Y. intercept. gen = Y l - (Slope. gen * X 2 );

/* Determine the Unknown Generator Point */

Generator. X[J] = (Y. intercept. perp - Y. intercept. gen) /
(Slope. gen - Slope. perp);
Generator. Y[J] = Slope. gen * Generator. X[J] + Y. intercept. gen;

for K=0; K<= Number.of.generator. points + 1; K = K + 1;

{
generate(G enerator. X[K], Generator. Y[K],
Generator. X[K+l],Generator.Y[K + l]);
}

} /* End Generate */

Figure 4.5. Detailed View of the Koch Algorithm (continued).

54
D. IMPLEMENTATION STRATEGIES
There are numerous ways to display the fractal shapes that the above
algorithm is capable of computing. The graphics primitives required are limited

to the standard initiation and termination commands coupled with the ability to

plot a point (or alternatively a line). Any raster graphics system, plotter or

similar technology suffices.

The algorithm can be extended to include:

- Online generator drawing to compose a generator relationship visually.

- Rotation in 3 dimensions (if your system has this capability).

- Variation of the inductive application of the generator by the inclusion of

randomness with respect to the generator constants.

Figures 4.6 through 4.9 represent a few of the shapes that this algorithm can

compute. Each figure has the generator data used to compute the shape and a

progression that shows the first two recursive iterations.

E. SUMMARY
The line intersection algorithm as it stands is not very useful for the

production of graphics images of realistically textured terrain. Its importance

results from its encapsulation of the essence of the non-random inductive fractal

method. This algorithm demonstrates the idea and intent of fractal functions and

their implementation within computer graphics. The potential fractal

programmer must throughly understand the salient parts of this chapter before

successfully attempting fractal images in three dimensions (i.e. before climbing the

mountain Chapter 5).

55
 t :

A Cloud-like Shape
Ini iat i ag Structure

Equi lateral Triangle

1st Iter&ti

GENERATOR DATA
Number of points = 1

Pointl :

Angle = 0.48S2 ra.ds .

Ratio = 1 .4222 2nd Itera.ti<.

Figure 4.6. A Cloud-like Shape

56
 : :

Boxes Ad Infinitum
Initiating Structure
Equilateral Triangle

1st Tterati

GENERATOR DATA
Number of points = 3
Pointl
Angle = O.OOOO rads
Ratio = 1 OOOO
.

Point2:
Angl e = O. 7854 rads
Ratio = l.OOOO
Point3:
Angle = 0.4636 rads
Ratio = 200.00

Figure 4.7. Boxes Ad Infinitum

57
 : .

An Exaggerated Koch Curve

Initiating Structure:
Equi lateral Triangle
GENERATOR/) I

INIT.

.w \

1st Iteration

GENERATOR DATA
Number of points = 3
Poiatl
Angle = O. oooo rads
Ratio = O. 6687
Point2:
Angle = o 7854 rads 2nd Iteration
Ratio = 1 OOOO
Point3:
Angle = o. OOOO rads
Ratio = 1 50O0

Figure 4.8. An Exaggerated Koch Curve

58
 .

A Plane Filling Curve

This plane filling curve requires
a slight modi fication to the dir-
ection of the generator app 1 i cm-
tion. This reversal is explained
on the next page

INIT. GENERATOR

GENERATOR DATA
Number o f points = 8
Point 1 :

Angl e = 1 OS3S rads

.

Ratio = 0.2000
Point2:
Angle = 1 OS3S rads
.

Ratio = 0.5000
Point3:
Angle = O. 7092 rads
Ratio = 2.0000
Point4:
Angle = O 3289 rads
.

Ratio = 4.0000
PointS:
Angle = O.OOOO rads
Ratio = O.SOOO
PointO:
Angle = O.OOOO rads
Ratio = 2.0000

-continued-

Figure 4.9. A Plane Filling Curve.

59
 Orientation

r~\ pi~^
of points
pa
r\
After Second Iterat:

/ ^S Angles = neg V^ \

rientation
Or \ / Orientation
of points \ / of points
Pi

Figure 4.9. A Plane Filling Curve (continued)

60
 V. FRACTAL GEOMETRY FOR GRA P HICS TERRAIN

One of the most widely recognized fractal images found in the literature is of

the mountain scene. This type of terrain modeling is perfectly attuned to the

fractal technique. The reason for this is that mountains are highly irregular
shapes, with a rough but consistent texture when viewed from a distant vantage

point. It is appropriate then, to introduce graphics terrain simulation techniques

through this model.

This chapter describes the theory and techniques of simulating mountainous

terrain with computer graphics. It provides the blueprint for fractal graphics

programming within R3 by providing general tools and a methodology that is

easily adapted to many other modeling needs.

A. MODELING MOUNTAINOUS TERRAIN

For the programmer who fully understands the essence of the method of

fractal programming introduced in Chapter 4, the movement into programming

in R3 is not difficult. The primary differences lie in the quantum jump in

computing resources that are required and the requirement to perform the

generator geometry in R3 versus R 2

. The theory and technique of fractals does

not change substantially.

Chapter 3 provided a rough framework to begin the coalescence of fractal

programming into a workable technique. We need to develop a number of tools

from that chapter and use standard compiiter graphics techniques to manage

those tools. To this framework, we add new fractal functions which provide the

texture of realism for our simulated mountain.

1. The Artist's Model

One way for an artist to build a physical relief model of a mountain is to

use a framework to provide structure to the model and a texturizing clay to

provide realism. The artist might use chicken wire on top of *mall boxes as the

frame with modeling clay as the texturizing element. His choice of clay is

predicated by the type of look that he wants to achieve. The chicken wire

CI
provides an inexpensive and disguised method to quickly build the mountainous
shape and structure. This method minimizes the cost and time to build up the

clay.

The artist continues the modeling process after the development of the

basic mountain shape to achieve hues and contrast in the coloration. He might
achieve this by the use of natural lighting to cast shadows or by a careful

painting of prominent features.

2. The Fractal Programmer's Model

There is very little in science that is truly new or innovative. We borrow

the essence of the above idea to guide us in developing a model for the discrete

computation of our two dimensional picture of the mountain. This section

describes the process intuitively and leaves the implementation details to later

sections.

a. The Lattice Control Structure

The pixel space that we developed in previous chapters can be

divided into discrete cubic units by use of a concept from mathematics called a

lattice (in our case we can view it as three dimensional graph paper). This lattice

serves as our controlling structure, the equivalent of the chicken wire structure

above. It is beneficial to build the lattice as a structure with well-formed

relationships, where the number of lines evenly divides the boundaries of the pixel

space and each line is a constant distance from its neighbors. By this method, we

do not have to store the lattice but can express it as a mathematical function of

the pixel space.

The lattice can be very useful in developing a rough approximation

of the mountain that we wish to model. This can be done in many different ways
but should result in a stick frame model of the mountain (a connected polygon

mesh like that of Figure 5.1).

The frame can be developed through an online graphics interface

that allows the programmer to select a ground level plane of the lattice and

provide a means to visually select points for the rough outline (essentially draw

the framework). This approach is useful when a particular shape is desired.

62
 A frame can also be developed using fractal functions to pervert the

lattice into a controlled random shape from a given plane of the lattice. This is a

powerful method that can be controlled via bounds on the random tools,

heuristics or discrete functional bounding of the fractal function. This approach is

most useful when a class of mountain shapes are required but no particular
mountain needs to be modeled, i.e. when random landscapes suffice.

Alternatively, the stick frame of the mountain can be determined via manual
(hand computation) means. This approach is tedious and limited in its flexibility,

and is not recommended.

b. Surface Texture via Fractal Functions

The next step in the creation of a mountain is to provide the

graphics clay to cover our stick frame model. This clay is a fractal function which

closes the polygons of the frame model with an inductive process that provides a
continuous pixel surface for the entire structure of the mountain.

The initiator/generator paradigm is used. The initial set of initiators

is the frame described above and the generator is a similar geometrical shape that

reduces in size continually until it becomes the size of a pixel and is mapped.
After the stick frame of the mountain is developed, this texturizing of

the surface becomes an automatic process that terminates when each geometrical
shape that makes up the framework is reduced to a continuous set of pixels in the

pixel set. At this point, the mountain exists in the pixel set (memory) but must
be provided color and light to bring it to life.

c. Hidden Surface Elimination

The pixel set has the entire stnicture of the mountain in memory,
but we can project only a two dimensional image of one plane onto the screen.

The fractal function which texturizes the surface does not concern itself with

local computations so many overlapping pixels are mapped to the pixel set.

There are two reasons then, why we need hidden surface removal (in this case

better referred to as hidden pixel elimination). First we have to eliminate the

back or hidden sides of the mountain by projecting only those pixels which are

visible along the axis of sight to the perpendicular planar surface of the display

screen. The second kind of hidden pixel removal is caused by mapped pixels

63
which were covered up by other recursive fractal descents either before or after

the pixel was mapped.

The removal of hidden pixels is greatly facilitated by the use of the

concept of the pixel set. In standard computer graphics hidden surface

elimination, the programmer is confronted with graphics primitives which are

functionally continuous Euclidean shapes. To effectively remove the hidden parts

of these shapes is in general very tedious and mathematically complicated. Since

the graphics programmer is shielded from the primitive -* pixel mapping, he is

26
functionally denied access past the simplifying abstraction of graphics Euclidean

primitives. The fractal programmer must have access to this level of the graphics

mapping and thus can use simple techniques to determine if a pixel is hidden or

visible.

The simplest and most economical means available to provide hidden

pixel elimination is through the use of Z-buffer algorithms. With this method, the

determination of whether a pixel is hidden can be appended to the pixel set

mapping process. The Z coordinate of a pixel that is to be mapped is checked

against the Z coordinate of the pixel currently in the pixel set at the same row
and column of the three dimensional array used to store the pixel set. If the pixel

is closer to the planar surface of the display screen, then the Z coordinate is

changed to reflect the position of the newly mapped pixel.

The Z-buffer approach, while powerful, does limit the fractal

programmer's flexibility. The axis of sight toward the mountain must be

determined prior to the fractal recursive process so that the determination of the

line through the (now two dimensional) pixel set is known. The Z-buffer becomes
an adjacency matrix to the pixel set and can retain information about forwardly
displayed pixels only. All information is lost about other pixels that were

computed in the fractal process. If another view of the mountain is required then

the entire pixel set has to be recomputed with a new axis of sight. If the fractal

function uses (non-tabular) random techniques then the mountain varies with

each view.

* The abstraction provided by Euclidean primitives is a powerful one when the alternative of
pixel mapping is considered. Without some powerful mapping tool (such as fractal functions), the
pixel level modeling process is in general very difficult.

64
 Most fractal pictures consume such vast computing resources that
only one view is computed for a given picture. As more requirements for graphics

terrain are determined, a more powerful method has to be used to retain all of

the computed pixels in the three dimensional pixel set. This method requires

that all pixels be stored in the three dimensional array previously described. The
hidden surface calculation can then be performed during the pixel mapping

operation or as a separate calculation that is performed after all fractal recursion

has terminated.

As specified in Chapter 3, the full array approach requires large

amounts of memory. This method, however, allows the computed fractal

mountain to become an entity that can be manipulated versus an instance of the

fractal mountain as above.

Both methods are viable but the latter approach provides more

flexibility for the programmer whereas the first approach is a response to the

economies of scale of data processing. As new architecttires are developed"' with

capacities geared toward fractal image computation, the first method can be
eliminated.

d. Illuminating the Mountain

If you stood on the dark side of the moon without illumination, the

mountains and craters of the moon would not be visible. They still exist however,

just as our imaginary mountain exists in memory. In order to visualize them, we

must illuminate them.
Illumination in computer graphics is achieved by varying the light

intensities of pixels displayed on the screen. The color mixture of these discrete

points determines the lighting effect that a viewer perceives. This perception is

not reality but another deception caused by scale and composition. A lighting

model then, is one which is able to abstract the essence of color from a real world

27
An Ideal architecture is one with a large main memory and parallel processing capabilities
with K processing elements (where K is greater than the maximum recursive descent distance).

65
object and transform that essence into a set of color values (intensities) that

accurately deceive the human eye via the graphics medium.

The literature on computer graphics contains many lighting models

with diverse approaches to the same problem. Many of these models (like those of

hidden surface) concern themselves with illumination of continuous Euclidean

28
surfaces and as such, are not directly germane to our study .

An object in space is a composition of basic elements. These elements

interact with the physics of light reflection to create the spectrum of light that

our eyes decode. In a graphics image, this process has to be simulated with

discrete lighting intensity values for each pixel. Thus, the illumination of the

mountain is a two step process; the fractal entities that are mapped to the pixel

set must be provided with a basic color, and these colors have to be highlighted

and dimmed by the lighting algorithm.

The basic color can be determined during the pixel mapping event of

the fractal recursion process or as a separate process prior to or in conjunction

with the lighting algorithm. This color can add realism to the picture through

heuristics which the programmer defines. Most mountains are composed of

different types of rocks and flora and these elements change at different altitudes.

This type of heuristic combined with some random control structure (i.e. to vary

the snow peak) can provide for improved realism (versus making the whole
mountain brown). The process of determining the basic color must be

accomplished prior to applying the lighting algorithm since the lighting algorithm

can only vary the intensities of an existing color 29 . Developing the process of

basic color determination is best accomplished through trial and error. It is the

artistic aspect of developing fractal mountains.

The general process of computer graphics illumination conr-rns itself

with casting shadows from one object to another given a direction from an

imaginary light source and with highlighting surfaces which are directly exposed

to the source. A surface is highlighted relative to the angle at which the light

"8
The Gouraud model (intensity interpolation shading) for instance

Since
Sine* so many diverse color models exist, the details of color representation are not
covered here

66
source's rays strike the surface. This poses special problems for fractal surfaces

due to their discontinuity at every point.

The process of illuminating a fractal surface is best aided by

divorcing the lighting process from the fractal computation process (except as

noted above). It is beneficial to view the pixel set as a collection of pebbles which
have size and position. This abstraction allows us to view the pixel as a

continuous space that can block light (cast shadows) and for which an angle of
illumination can be determined (usually in conjunction with neighboring pixels).

A well formed fractal mountain surface is completely connected (no

space between adjacent pixels in the pixel set). Thus the surface can also be

viewed as a continuous (while very rough) surface where reflected light can be

cast from or to adjacent pixels.

One lighting model which fits the fractal process is the Torrance-

Sparrow model [Ref. 8:pp. 578-579].

This model views an object as a collection of facets which is each a perfect reflec-
tor (i.e. does not absorb light). The orientation of each facet is given by the
Gaussian probability distribution function (i.e. the smooth surface of the Eu-
clidean object is roughed by the Gaussian relationship). The geometry of the
facets and the direction of light (assumed to be from an infinitely distant source.
so all rays are parallel^ determines the intensity and direction of specular reflec-
tion as a function of the light source intensity, the normal to the average sur-
face, the direction to the light source and the direction to the viewpoint.

This model has to be modified to adapt to the fractal set method. In the fractal

method, there is no need to rough the surface to provide reflection because the

surface is by design roughly textured. A method of assigning planar front.9 to

each pixel space has to be determined and the geometry of connecting these

fronts identified. With these modifications to the lighting model, each individual

pixel's color intensity can be modified for the increase in intensity associated with

the light which falls upon it.

The model also allows diffuse reflection (light reflected from one

object to another) which is critical to bring out clarity of the fractal image. For

further information on the model the reader is referred to the reference.

67
 e. Summary of the Fractal Mountain Paradigm
To summarize the methodology we can view the process as a five

step process:

- Build the initiator framework or stick frame model of the mountain.

- Give the frame's surface texture with fractal functions.

- Remove hidden surfaces (pixels) from the display.

- Illuminate the surface with lighting algorithms.
- Project the surface to the screen.

B. FRACTAL TOOLS FOR TERRAIN MODELING

The tools presented in this section can be used in the creation of fractal

images within R3 . The list provides a basic set of programming tools to guide the

creation process.

1. Equations of the Lattice

The lattice (or controlling structure) can be very useful to the graphics

programmer to implement heuristics or bounding functions on the essentially

random progression of the fractal figure. The graphics programmer may wish to

limit the growth of the mountain by implementing a ground level plane of the

lattice and a maximum height that the mountain can obtain. He accomplishes
this by arbitrarily assigning another plane of the lattice as the upper bounding

plane. The height of the mountain can then be checked during any level of the

fractal recursive descent against this fixed plane. The programmer can then clip

the height by adjusting the random equation that controls the upward trend to

tend towards the ground again. This is an example of a heuristic applied to the

fractal recursion that controls the external qualities of the function.

A fractal programmer can use the lattice to assign the initial colors to

the mountain via a user designed set of rules. The lattice aids the user in the

implementation of the rules by giving reference points for inclusion of branching

conditions (tree line to snow line etc.) and can be used in conjunction with

decision weights to add a varied transition between textures. The determination

of the initial color of a mapped pixel is usually a controlled random process, one

of the primary methods of control being the lattice or some derivative thereof. As

68
the height of the mountain increases (lattice level), it becomes increasingly more
likely that it will transition to another texture. This can be controlled by adding

the lattice level as a factor to the rule that decides color.

An example of a potential lattice equation and how it might relate to the

pixel space is demonstrated in Figure 5.1. The actual lattice has been extended

from the pixel space in order to visually demonstrate how it relates to the pixel

set. In actuality, this is not the case. The lattice coincides with the boundaries of

the pixel space. Although the lattice can have a one-to-one relationship with the

pixel space, this defeats the purpose of the lattice (macro control). By grouping
cubic sets of pixels into a well-formed relationship, we can better implement
heuristics and bounding functions.
As a lattice example, consider a pixel space that is created by abstracting

the real world coordinate space for our mountain as described below. We desire a

real world space to be a cubic area established by the box 20.000 ft.

(x coordinate ) by 15,000 ft. (y ) by 20,000 ft. (z ) This can be sectioned into a

lattice by establishing the increment of distance between adjacent lattice points

30
to be 1000 ft. and establishing the corner lattice point as (0.0, 0) .

The mapping function between the lattice and pixel space is then

straightforward. The size of the pixel (recall equation from Chapter 3) is

8
: -=20 ft. and the lattice cubic sections contain 10 cubic feet or
1,000 ft.

equivalently 75,000 cubic pixels.

The ground level can then be identified as the 2000 foot level and the

bounding height can be assigned a level of 10500 feet. If we wish, we can make
the bounding heuristic more realistic by sectioning the lattice into mountainous

areas, each having different bounding levels.

A completely arbitrary set of dimensions, increments and points.

69
The Lattice Control Structure

LATTICE
FRACTAL MOUNTAIN STRUCTURE
PIXEL SET

Fig 5.1. The Lattice Control Structure

70
 2. A Fractal Function for Contouring Mountains

The usual method for contouring mountains uses a randomized variation

of the mid-point displacement method introduced in Chapter 4. The planar

structure is typically the triangle

31
imbedded in R3 . The basic methodology is

demonstrated in Figure 5.2. Figure 5. 2. a shows a triangle with its first iteration of

mid-point displacement. This process continues until all triangles have reached

the desired level of precision. One completed structure is demonstrated in Figure

5.2.b. The precision is typically lower (pixel level) than that demonstrated in

Figure 5.2.b but it was terminated at a higher level to better demonstrate the

idea. Random techniques (described below) are used to produce the relatively

accurate picture of a mountain frame as depicted in Figure 5.2.c.

In practice, the random techniques are implemented with the mid-point

displacement function during the fractal recursive descent. The random
techniques provide local disorder to the fractal function which provides the

computational structure. Results have shown that very little randomness needs to

be applied to the regular structure of Figure 5.2.b to achieve satisfactory results.

The mountains created for the film Star Trek: The Search For Spock used a

limited random number look-up table consisting of fewer than 300 entries [Ref.

-]

3. The Geometry for Mid-Point Displacement

The general approach to building a fractal shape as illustrated in Figure

5.2.c is to use the algorithm of midpoint triangle displacement combined with a

randomized displacement along the normal to the X-Z plane of a cartesian three

space coordinate system. A recursive procedure which computes this relationship

requires as inputs the points of the triangle. It computes the midpoints of each
line of the triangle and inscribes a triangle inside of the initiating triangle by

connecting each midpoint, Figure 5. 3. a. This process yields four triangles

coincident with the plane of the initiating triangle. When we fix the X-Z normal

at any of the midpoints, we can displace the midpoint by a discrete distance

1
Any regular structure suffices; the triangle is easy to use and yields very satisfactory
results.

71
 .

The Triangular Midpoint

Displacement Technique

Determine Connect Connect Connect

Midpoints Ml to M2 M2 to M3 M3 to Ml

Fig 5. 2. a The 1st Iteration of Midpoint Displacement.

The Random Structure was rotated

-30 degrees around the X axis to
accentuate its texture

Fig 5.2.b Completed Structure. Fig 5.2.c Randomized Version

Figure 5.2. The Triangular Midpoint Displacement Technique

72
along the normal and determine a point, Figure 5.3.b. Since the normal is to the

X-Z plane, it is sufficient to simply modify the Y coordinate according to a

positive or negative value. This is equivalent to displacing the midpoint along

the X-Z normal up or down. We perform this displacement to each midpoint

normal and replace the midpoint with these new points. This yields a new
structure that still consists of four triangles but with each coincident with a

different plane, Figure 5.3.c.

a. Midpoint of a Line in R 3

The determination of the midpoints of the lines of the initiating

triangle is a simple process that uses the equation of Chapter 4. and fixes the

generator ratio constant at 1. This simplifies the general equation of:

Xj + (Generator. ratio. constant x X 2 )

1 + Generator. ratio. constant

to the well-known midpoint relationships of:

X, + x 2
Xmid
2

Y, + Y
Y mid
*
2

z, + z2
Z m id

The above equations completely determine the midpoints of the lines formed by

each endpoint of the initiating triangle.

b. Displacement along the X— Z Normal

The process of displacing the midpoint along the X-Z normal is a

simple one. We need a factor such that the displacement can obtain a varied

magnitude. This is best aided by the inclusion of a random variable as a multiple

of some scaling factor that is added to the Y coordinate of each computed

midpoint. This process is demonstrated in the following code segment:

Randvar = getrand(Seed);
Point l[y] = Pointl[y] + (Scale * Randvar);

73
Four Triangles are
created for each
triangle initiator

Figure 5. 3. a. Triangular Midpoint Displacement

Computed displacement

After 1st displacement

X-Z plane

Figure 5.3.b Displacement Along the X-Z Normal

Original generating triangle

Negative Y displacement

Positive Y displacement

Figure 5.3.c Completed Random Fractal Triangle.

Figure 5.3. The Random Midpoint Displacement Technique

74
 A valid question is, why the normal to the X-Z plane? There are

three good answers to this question. Using the normal to a fixed plane simplifies

the computation (eliminates the need to perform planar computations at each

recursive division). It also is generally the direction that we want the mountain

to grow. The most important reason however, is related to the gapping problem

(described below). With a fixed direction for displacement, there is no need to

communicate the direction of displacement along the normal between adjacent

side computations. The recursive levels that compute adjacent sides are

functionally discordant. It is demonstrated below that the solution to the gapping

problem (inconsistent random numbers) which creates the need to communicate

along discordant recursive levels is algorithmically difficult to solve and thus

should be avoided.

c. The Gapping Problem

One problem exists for the midpoint displacement procedure which

utilizes a random displacement along the X-Z normal line. It is indirectly caused

by the data locality aspect of the inductive process of the recursive fractal

descent. The problem exists when two adjacent sides of two adjacent triangles are

not displaced with the same value. Each side is computed during independent
levels of the recursive descent so there is no practical method to communicate the
random numbers for the displacement.

The gapping problem is illustrated in Figure 5.4. For the two

triangles that are extrapolated from the structure, there is an unknown

relationship that is the refndom variable used to displace the common midpoint.

If triangle A uses Rand = 0.3 and triangle B uses Rand = -1.07, then the

displacement for each adjacent midpoint (which are at the start coincident) is

skewed in the opposite direction. This creates a gap in the fractal landscape that

will (in all likelyhood) not be filled by other fractal shapes from neighboring

triangles. We need an algorithm which can insure that each midpoint (which is

always shared by two triangles) has the same displacement along the normal to

the plane.

75
 d. Solving the Gapping Problem via Random Tables

The solution to the gapping problem is straightforward if the

programmer adopts the random number table as his random function

implementation. The goal is to match adjacent triangles with a seed or

displacement within the random table so that the random number returned is

equivalent for each coincident midpoint.

There exists a symmetry within the triangle of Figure 5. 5. a that

allows such an approach. Ideally we want the point M a to be displaced by the

same magnitude when triangles T 2

and T4 (highlighted by textures) compute

their random numbers for M a . This can be facilitated by the inclusion of a table

seed for each recursive call to the midpoint displacement routine and by rotating

the orientation of the midpoint triangle (the triangle created by the three
computed midpoints) labeled T4 in Figure 5. 5. a. This rotation is performed in

relation to the random table and not in relation to the Cartesian space. It is

accomplished by adjusting the order of the points in the recursive call.

The order of the points for triangles Tj T2 and T3 are as described

in Figure 5.5.b and for T4 as described in Figure 5.5.c. All four triangles generate

a recursive sequence and use the same seed to the random number table. The
random numbers retrieved from the table must observe the order of assignment
that is demonstrated in Figure 5.5.d. For example, the line segment formed by

the first two points (Pj and P2 ) input to the midpoint displacement routine

determine the midpoint Rj. This midpoint is assigned the random displacement
from the table corresponding to the entry seed. The next midpoint retrieves the

table entry corresponding to seed + 1 and so on.

If this technique is followed the sequence of random numbers will

match-up as demonstrated in Figure 5.5.e. The recursive calls correspond to the

code segment in Figure 5.6.

76
 The Gapping Proble m
Gene rated random structure

Ini tinting Triangle

ADJACENT TRIANGLES SHARING

A COMMON MIDPOINT

common midpo i nt

AFTER MIDPOINT DISPLACEMENT

Gap created by discontinuous

random displ acements along
the Y axi s
.

Shared common midpoint

Figure 5.4 The Gapping Problem

77
 b da
e

Solving the Gapping Problem

t3

Shared ton midpoint

mm
/ Tl >
Ip^iik
a \
Figure 5 5 . .

Figure 5 . 5 .

Given that the input to the random table is seed=i

Figure 5 5 . .

i+l

Figure 5 5 . .

Figure 5.5. A Solution to the Gapping Problem

78
Main()
{
LOAD THE INITIATING TRIANGLE
Seed = 1;

frac_triangle(P j ,P 2 ,P 3 ,Seed)

frac_triangle(Pj,P 7 ,P 3 ,Seed)
{

DETERMINE DISTANCE BETWEEN ENDPOINTS OF AN INITIATOR

If (DIST < Pixel. length)

Plotj>ixel();
return;
}

COMPUTE THE MIDPOINTS (M .M .M 3 1 2 )

ADJUST THE Y COORDINATE FOR Mj Using Randtable(Seed)

ADJUST THE Y COORDINATE FOR M 2 Using Randtable(Seed+l)
ADJUST THE Y COORDINATE FOR M 3 Using Randtable(Seed + 2)

Seed = Seed + 3;

/* Triangle T, */
frac_triangle(M .P 2 ,M 2 ,Seed)
1

/* Triangle T2 */
frac_triangle(M 3 .M 2 ,P 3 ,Seed)
/* Triangle T3 */
frac_triangle(P .M, j .M 3 ,Seed)
/* Triangle T4 '/
frac_triangle(M 2 ,M3 ,M 1 ,Seed)
}

Figure 5.6. An Algorithm for the Midpoint Displacement Technique.

79
 4. Random (Stocastic) Fractals

One common complaint about computer graphics images and animations

is the artificial perfection of the displayed shapes. Our mind subconsciously rebels

against the order that is displayed, our expectations about the rough reality of

nature are not satisfied. The use of randomness in generating fractal images is

necessary to approximate the observed disorder of nature. An example is the

Koch curve. Although it resembles a snowflake, it lacks the realistic look that

experience trains our eyes to see. In a mathematical sense, the Koch curve is

beautiful; as an approximate to nature it lacks appeal.

To approximate the rough texture of nature, we are forced to modify the

well-behaved mathematical relationship of the initiator^generator in a controlled

manner to add variety to our computed image. This modification is usually by

the inclusion of a random variable into the control structure within the fractal

equation. The random variable must exhibit restraint. It cannot be allowed to

vary wildly without structure.

One of the most appealing random functions which provides very

satisfactory results in fractal images is the normal distribution 32 . The normal

distribution (as opposed to a uniform distribution) approximates the expected

local disorder in nature (at least experimentally).

a. The Normal (Gaussian) Distribution

The normal distribution is used throughout the natural sciences for

many applications. It was first derived as an empirical result of the observed

error about a true value that normally occurs when measurements are taken of a

natural event. The symmetry that was observed from error measurement and
sampling suggested that there was a natural order to such observations. These

empirical results spurred natural scientists and mathematicians to try to fit a

curve to the observed graph that behaves as probability requires (i.e. the sum of

the area under the curve equals unity). Many of the early scientists

32
Often referred to as the Gaussian distribution, the normal distribution is the standard bell
curve to which every student is accustomed.

80
referred to the normal distribution as the law of error in deference to its roots in

experimental natural science.

Many functional characterizations of the normal distribution were

developed 33 but credit

, is usually attributed to Carl Frederic Gauss [Ref. 9:pp. 1-

11] who formulated a least squares approach, published in 1809 in Theoria Motus

Corporum Coelestium. The form of the normal distribution was not finalized

until the early 20 century.

We take our definition of the normal distribution from [Ref. 9:pp 18].

refer to Figure 5.7.

Definition:

The probability density function of a normal random variable X is given by:

ft is the mean, a is the standard deviation
and a 2 is the variance.

f(x: M V "Ik
•exp
2a
2

where -ao< x < oq -oo< fi <oq and ct>0

The Normal Distribution

SCALE:

horz vert

Figure 5 7. The Normal Distribution

ss
Many mathematicians can lay claim to founding the normal distribution, most notably: Pi-

erre Simon de Laplace and Abraham de Moivre Ref. 9:pp. 11

81
 This definition is the general case of the normal distribution. We are

interested in the behavior of the function and need a practical way to determine

a random number that we can use in the parameter of the normal to the plane in

the geometrical relationship described above. To facilitate this, we simplify the

general normal distribution to the well known standard normal distribution,

illustrated in Figure 5.8. The standard normal distribution function is the special

case where /i = and a = 1. This reduces the general equation to the simplified

equation:

i x"
f.M = exp
'2* T
The above functions describe the behavior of a normal random

variable. We need a function that returns values from that function which will

observe the period of the normal distribution. This means we need a string of

real numbers over an assigned range about a mean that will observe the

frequency of the normal distribution.

The Standard
^Normal (Gaussian)
SCALE: Distribution
Wz a 5
vert
/*= 0.0
(7= 1.0

N. (x) = 7±= e 2

>/iT

99.37%-

Figure 5.8. The Standard Normal Distribution

82
 b. Standard Computer Random Functions

Some computer systems provide a random number generating

function which observes the normal distribution. If this is provided, then it can

be used directly (after scaling) as a parameter to displace the Y coordinate in the

geometry of the normal to the midpoint displacement as described above.

Many computer systems only provide a random number generating

function which is uniformly distributed over an interval of integers. This is a

pseudo- random number. Such a function, when given a seed, will produce a

sequence of numbers distributed over the fixed interval defined by that system.

The interval is typically proportional to the maximum integer defined in the

compilers of the system. A normal distribution routine must then be defined that
transforms the uniform random numbers into random numbers which behave

according to the standard normal distribution function.

There exist transformation functions that take a uniform random

variable distributed over the interval [0,1] into an approximate normal random

variable over — oo < x < oo . This requires the uniform random variable to be

mapped into the interval [0,1] and then transformed by the normal

approximating function.

To transform a uniform random variable distributed over an interval

[0,max_int] while maintaining the distribution density, requires the following step:

UNF- 0ma ^ n
UNF I0>1]
=
max int
(

One commonly used function that transforms uniform random

variables into normal random variables is found in [Ref. 9:pp 49]. This function

uses two uniform variables from [0,1], denoted UNF, and UNF 2
. and computes

two normal random variables, denoted NORMj and NORM,.

NORM, -
>
/(-2log eUNF 1 cos(2*rUNF 2 )

NORMo = v /(-2log e UNF sin(2?rUNF 2 ) 1

S4
This is how most standard system provided computer subroutines perform the operation

83
Appendix B contains a C UNIX routine that implements an algorithm to

compute the um'/orm[0,l] -» normal [— oc,+oo] transformation.

A programmer must be very careful when dealing with random

number generators from standard system subroutines. These routines vary widely

and can provide good to barely adequate results. When the normal

transformation routine is written, the programmer must verify experimentally

that his function adequately models the normal distribution. This process is

illustrated by Figure 5.9. Appendix B also contains experimental results which

verify the transformation.

The purist may not accept the results displayed in Figure 5.9 as an

accurate transformation (there appears to be a skew to the negative direction).

We must remind ourselves that we are trying to approximate the roughness of

nature and minor random skewness will not deter us. If the programmer demands
a better approximation, it is a simple process to expand the sample space of the

test and build a table with exact proportions by selective deletion of skew

density.

c. Random Functions versus Table Driven Methods

The application of a random modifier in the midpoint displacement

technique can be achieved via two methods.

- By invoking the above function iteratively as a variable.

- Or by a variable returned from a table lookup operation from a random
table.

The choice of which method to use depends on the programmer's application but

each has its ramifications.

In general, the table lookup operation is considerably faster than the

functional method but must by its definition limit the amount of randomness it

contains. The major issue however is the need to reproduce a figure under some
requirement for fixed terrain. This issue was the driving force for Loren Carpenter

from Lucas Film in determining that he needed to use a table driven method to

produce the planet images for the film Star Trek: The Search for Spock [Ref. 7].

He had to be able to fix a space where the images of the actors could be imposed

84
 : . :

Experimental Results for the

Computer Generated Normal
Distribution
Sample Space:
500 Random Events 115

Experimental Reiulti
With Normalised
a.

Uniform Distribution
Transformation Equations Over [0.0,1.0]
From Uniform Distribution
Over [0,1] to the Normal 87
so-
Distribution Over -oo<x<oo
40- issgj 8
30-
N2 =y (-2lnU1 ) * cos(2KU
2)
2CT
::•
:•:

1Q-
N 2 = / (-2lnUJ ) * sin(2tU ) 59
58
o i o
Samp 1 e Spac e
SOO Random Events

29

.21 i .

18
12T
T' T T

' ,i ,

ri
i

r"i "i i i"i"

-3.0 O.O 3.0
Figure 5.9. Experimental Results for a Computer
generated normal distribution over [0,1]

85
onto the fractal images and could not allow the fixed space to change with each

frame computed. This is the major advantage of the table method. By retaining

a seed to a table of random numbers, you can reproduce the sequence of

displacements along the normal during the fractal recursive descent.

When you consider the existential qualities of randomness you are

confronted with basic questions about determinism and order in the universe. It

is not at all clear which rules chance. In any case, we can deceive perception

with a relatively small table of random numbers.

The question of how much randomness is enough to provide for a

visually appealing texture is not completely clear. In [Ref. 7] Smith demonstrates

a variety of shapes computed with the same algorithm of Figure 5.6 using

random number tables of different sizes. He demonstrated that as few as five

numbers can suffice to provide enough local disorder to give the viewer the

acceptable texture of a mountain. If the mountain segments are viewed at the

correct perspective and scale, this perception is clearly felt. A trained

mathematician would find the five element mountain statistically unappealing

however. A true stochastic construction requires a continuous random function

rather than a discrete table method. As long as the goal of our computations is

merely to deceive the graphics viewer, it suffices to use the random number table.

The table must be large enough to provide for an appealing textural perception.

A complete C program that computes a triangular mountain segment using the

random displacement midpoint technique is contained in appendix C.

86
 VI. SHORT CUTS TO MOUNTAIN SHAPES

Since the fractal mountain computation (the full approach with hidden

surfaces etc.) is so costly in terms of resources, it is important for us to consider

shortcuts that can lessen this burden. This is best realized by utilizing the hidden

surface and curve fitting capabilities that are provided on some advanced

graphics systems.

Our goal is to match the well known bicubic surface procedures with the

structure computed by the simple fractal algorithms. This is best accomplished

by modifying the triangular midpoint displacement technique and using a

35
rectangle as the basic geometric building block. Most of the cubic surface

algorithms use the rectangular structure as their basis, so it is easier to adapt

them to our fractal structure.

When the fractal algorithm of Figure 6.2 has its computations terminated

before reaching the level of pixel size, it yields a connected rectangle structure

like the one shown in Figure 6.3. This structure is a connected Euclidean

structure that can be used as a base on which other algorithms can be applied.

Cubic equations can fill the polygons to an arbitrary precision and standard

hidden surface algorithms can eliminate the hidden sides of the computed

surfaces Simple lighting algorithms can be applied to the computed surface to

36
achieve a realistic lighting effect . This is how Voss and Carpenter created their

fractal surfaces in [Ref. 7].

A. RECTANGULAR MIDPOINT TECHNIQUE

Modifying the triangular midpoint algorithm of chapter 5 is a straightforward

process that introduces no new mathematics or difficulties. It consists of a

procedure to split the midpoints of each side of the rectangle and a procedure to

find the center of the rectangle. From these five points, we construct four scaled

55
We actually use a non-planar four sided polygon. We refer to the basic structure as a rec-
tangle to simplify the terminology.
6
Couraud shading for example.

87
rectangles, as demonstrated in Figure 6.1. The five shared midpoints of the

generated rectangles are then displaced along the normal to the X-Z plane

according to a random Gaussian value. This process is exactly the same as for the

triangular algorithm of chapter 5. The gapping problem still exists and this

requires an algorithm to rotate the rectangle relative to the random number table

and the starting seed to insure that adjacent midpoints are displaced relative to

the same random number. The basic methodology is displayed in Figure 6.1, the

algorithm is contained in Figure 6.2 with sample results in Figure 6.3.

88
 Computed
Four Rectangles Midpoints
are created for
each rectangu-
lar initiator.

Computed displacement

Y axis

'X-Z plane

Negative Y
displacement

Positive Y displacement Original generating rectangle

COMPLETED FRACTAL RECTANGLE

Figure 6.1. The Rectangular Midpoint Displacement

89
Main()
{
LOAD THE INITIATING RECTANGLE
Seed = 1;

frac_rectangle(P j ,P 2 ,P 3 ,P 4 ,Seed)
}

frac^ ectangle(P ,P 2 ,P3,P 4 ,Seed)

,

DETERMINE DISTANCE BETWEEN ENDPOINTS OF AN INITIATOR

If (DIST < Pixel.length)
{
Plot joint ();
return;
}

COMPUTE THE MIDPOINTS (M ,M ,M 3 .M ,M C 1 2 4 )

ADJUST THE Y COORDINATE FOR M, Using Randtable(Seed)

ADJUST THE Y COORDINATE FOR M 2
Using Randtable(Seed+l)
ADJUST THE Y COORDINATE FOR M Using Randtable(Seed+2)
3
ADJUST THE Y COORDINATE FOR M Using Randtable(Seed+3)
4
ADJUST THE Y COORDINATE FOR M c Using Randtable(Seed+4)
Seed = Seed + 5;

/* Rectangle R 2
*/
frac^;ectangle(P 1 ,M 1 ,M c ,M 4 ,Seed)
/* Rectangle R 2
*/
frac_rectangle(M 2 ,M c ,M ,P 2 ,Seed) 1

I* Rectangle R 3 */
frac_rectangle(M 2 ,M c ,M 3 ,P 3 ,Seed)
/* Rectangle R 4
*/
frac_rectangle(P 4 ,M 3 ,M c .M 4 ,Seed)
}

Figure 6.2. An Algorithm for the (Rect.) Midpoint Displacement Technique.

90
 Rectangular Mountain Fractal

Generating rectangle

The recursive termination event was 1/2 inch

Figure 6.3. An Example of the Rectangular Mountain Fractal

91
B. PARAMETRIC CUBIC SURFACES
A complete description of parametric cubic surfaces is too involved to be

described in this study. The theoretical basis of cubic curves is not directly

applicable to fractal geometry. For a complete description refer to [Ref. 8:pp.

514-536]. If the reader is already familiar with cubic curves and their derivations,

he can skip by the section on cubic curves to the section that details the

application of cubic surfaces. For any reader who has not been exposed to the

derivations of parametric equations which yields cubic curve computational

engines, it is recommended that he read the following section so that he may gain

insight into the mathematics of cubic surfaces. Detailed knowledge of cubic

curves is not a prerequisite to the successful use of cubic surface fitting engines

with respect to fractal surfaces. It is helpful, however, to understand the

underlying mathematics whenever canned equations are used.

1. Cubic Curves

The general method of cubic curves has as its basis that any continuous

curve in R can be expressed in parametric form. This form relates the points

x,y,z with a parameter t such that as t varies within some range of values the

equations solve for unique points on the curve. Specifying two endpoints and two

control points of a segment of the curve allows us to define certain constraints to

be applied to the parametric equations. These constraints allow us to manipulate

the parametric form of the equations to yield a simple vector product definition

of that segment. Once this vector product is established, we can solve for points

on the curve by picking discrete values of t and solving for x. y,z in turn. This

yields a discrete approximation of the curve that can be as precise as needed,

a. Parametric Cubic Equations of a Curve

A parametric cubic curve is one for which the points in R3 (x,y,z)

are each represented as a third-order (cubic) polynomial of some ;>arameter t.

Because we deal with a finite segment of a curve, we limit the range of the

t may vary between and 1 for example.

92
parameter t to the range, ^ t ^ 1. This yields the equations:

x(t) = a xt 3 + b x t 2 + c xt + dx
y(t) = a t3 + b t2 + c t + dy
y y y

z(t) = a z t 3 + b,t 2 + e,t + dz

Each equation can be expressed as a vector product as x(t) is below:

3
x(t) = [t t
2
t 1]

This vector product separates the distinct parameters of the parametric equation

into the unknown coefficients of x(t); [a x 'D x c x^x] an(^ ^e parameter t that we
wish to manipulate. Through this separation, we are able to manipulate them as

algebraic entities. If you multiply the vector product out. you find that the

vector product is equivalent to the parametric equation that precedes it. Denote

this product as x(t) = TC X where

3
T = ft t
2
t 1

and

The vector T is the same for x(t). y(t) and z(t).

We now establish constraints (as a set of control points) for the

equation x(t) evaluated at the bounds of the range of the parameter /. (i.e. t=0
and t = l). We consider four equations of x(t) and its first derivative x'(t) where

these boundary conditions yield four known points.

x(t)=TC x when ; evaluated at=0. -> x(0) =[000 l]C x

x(t) = TC x ; when evaluated at=l. -» x(l) = [l 1 1 l]Cx

93
and since the first derivative of x(t) is:

x'(t) = [3t
2
2t 1 0]C X = T'C X

x'(t) = T'Cx ; when evaluated at=0, -> x'(0) = [0 1 0]C X

x'(t) = T'Cx ; when evaluated at=l, - x'(l) = [3 2 1 0]C X

We now have four equations that can be grouped into a vector product:

x(0) 1

x(l) 1111
x'(0) 10
x'(l) 3 2 10

We recognize that x(0) and x(l) are the endpoints of the curve

segment and x'(0) and x'(l) are components of the tangent vector at the

endpoints (y'(t) andz'(t) are the other components). With this knowledge we

are able to solve the left hand side of the equation above. These points (that we
38
call Pj through P4 ) are the control points that we establish for curve fitting .

For a given curve segment the control points are fixed. We rewrite the equations

above with respect to these known control points:

1 1 1 1

3 2 1

Denote this equation as:

Gx = MC X
The matrix Gx is often referred to as the geometry of the cubic curve and M as

the basis.

This equation has the 4 by 1 row vector C x as the only unknown.

The elements of the Cx vector are the parameters (a x ,b x .c x .d x ) from the

If we establish ourselves as servers then these four points are the user's input to our rou-
tine.

94
parametric equations. We can solve this equation for these parameters and

establish the parametric equations with the only unknown being the parameter /.

The parameter t can be discretely varied over its range of ^ t ^1. providing a

set of points on the curve. It is through these constraints that the control points

control the parametric equations and produce an equation that can produce a

discretely sampleable curve segment in three space. Solving the equation for Cx is

straightforward:

CX =M 1
GX
39
Substituting Cx into the equation for x(t) yields :

x(t) - TM 'G x

Similar arguments yield the equations for y(t) and z(t):

y(t) - TM-'Gy
z(t) - TM'G,
The matrix M _1 is constant for all three equations and is usually

denoted by the type of surface that it relates to Bezier — Mb , Hermite — Mh

etc. It is through the control points and their interaction with the constraints

that the models Bezier, B-spline, Cardinal Spline, Ferguson (Hermite or Coon's)

surface etc. modify the parametric equations and provide different curve fitting

engines.

For each model, the matrix M mode j is constant throughout all

computations. To use the model requires the determination of the control points

(in conjunction with how they relate to the curve) and a vector multiplication
engine. Since vector pipeline computations are ideally suited to computers, this

method becomes a fast technology for curve fitting with an intuitive appeal for a

programmer.

We have just determined the Hermite model equation for x(t).

95
 b. An Example: Bezier Cubic Curves

We consider the model called Bezier [Ref. 8:pp. 514-536]. The Bezier

model defines the position of the curve's endpoints and uses two other points (not
on the curve) which define tangents at the curve's endpoints (by the line segment

joining the tangent points to the endpoints).

The matrix M is derived by setting the following constraints (see

Figure 6.4). One endpoint of the segment is located at P^

x(0)=P 1

The other endpoint is located at P4 :

x(l)=P 4

The line segment from Pj to P2 defines a tangent at P 1

such that x'(0) relates

to the points Pj,P 2 as below:

x'(0) = SfPr-PJ

And similarly for the tangent at P4 defined by P 3 ,P 4 :

x'(l) = 3(P 4 -P 3 )

Solving for Cx in terms of M b yields the cubic Bezier matrix as:

1 3 -3 1

3 -6 3
G,
3 3

Hence the equation for x(t) is:

1 3 -3 1
Pi

3 -6 3 P 2
s z
x(t) = [t t t 1]
3 3 P3
1 P.

96
 Bezier Curve
Tange nt, defined by
P P line segment.

Tange nt defined by
P P „ line segment.
4 <1

Figure 6.4. An Example of a Bezier Curve

The process of creating a Bezier curve given the above parametric

cubic engine is a simple process of computing discrete points on the curve by

substituting values along the range of t and fitting the curve by connecting each

point with a line segment. This provides an approximation to the curve that can

be processed at an arbitrary precision by incrementing St with smaller and

smaller lengths.

The process of shaping a curve is accomplished by increasing or

decreasing the two endpoint tangents formed by the four control points. It can be

viewed intuitively by thinking about each tangent as a force which pulls the

curve in the direction of the tangent until the force from the other endpoint

overcomes the original at the midpoint. The two endpoint tangents work against

one another proportional to the distance of St from each endpoint.

2. Bezier Surfaces

Extending the above method to cubic surface sections is accomplished by

adding a new parameter s that we vary from 0^5 ^ 1 as we did with the

97
parameter t in cubic curves. The connection between cubic curves and surfaces

can be made by fixing one parameter and varying the other over its range. This

yields a cubic curve. The equation is of the form x(s,t) and is written as:

x(s,t) = a n s 3 t 3 + a 12 s 3 t 2 + a 13 s 3 t + a ]4 s 3

T ^41' t ®42 "*" ®43 ' ®44

Written in the algebraic form:

x(s,t) = SC x T e

where S = [s
3
,s
2
,s,l], T = [t
3
,t
2
,t,l] and T fc
is the transpose of the matrix T.

The complete algebraic manipulation of the equation to arrive at the

equation below is similar to the curve process as described in the previous

section. Its details are covered in [Ref. 8:pp. 524-536]. The equation for a Bezier

surface patch is:

x(s,t) = SMbQ^^T*
where Mb is the same matrix as in the curve equation, Mb 1
is its transpose and

Qx is the x component of sixteen control points of a surface patch. Bezier

surfaces are intuitive in their appeal and serve the fractal rectangular mountain
well. To apply the technique to the mountain of Figure 6.3 requires the

application of a routine that takes the non-planer four sided shape of a computed
initiator and develops a connected sixteen point figure as illustrated in Figure 6.5.

The inclusion of a Bezier subroutine at the recursive termination event after this

figtire is developed matches the sixteen point figure with a smooth curve. To

achieve edge continuity requires that adjacent sides have the same four points in

proper juxtaposition in the sixteen point matrix. This is also demonstrated in

Figure 6.5. Bezier surfaces guarantee such continuity.

98
 — ,

Bezier Surface Patches

— —— — — —_—
- —
, _ ,
_ . , __—, ——
. . __ —

p ^***-* Bi m Shared edge points

Qx for R, =

Qfor R, =

Figure 6.5. Matching Bezier Surface Patches to the Fractal

Rectangular Mountain Structure.

99
 VII. CONCLUSIONS

A. DIRECTIONS FOR FURTHER STUDY

Fractal geometry as an area of research is very new. Because of this, there is

a great need for refinement and exploration. What is known needs to be refined

into a set of workable techniques with reasonable, simple terminology as its root.

The areas that are unknown need to be explored intrepidly. With this goal in

mind, the following paragraphs quickly review some areas of prospective research.

The reader is invited to explore their potential.

1. Development of New Fractal Functional Methods

The current tools of fractal functions are tentative and limited in their

ability to yield insight. New applications of the recursive initiator- generator

paradigm are waiting to be discovered. This area of research is especially good for

the graphics programmer since the graphics medium is currently the best method
for fractal experimentation. As these new functions are developed, they can be

shared, yielding a glossary of modeling functions that can be molded into a

cohesive theory 40 . Related to this is the need to develop a functional language

(within the language of mathematics) of fractal geometry to aid in the

communication of ideas and in the eventual coalescence of the theory.

2. Fractal Lighting Model

The current state of the art in computer graphics lighting models lacks a

complete model for the pixel set paradigm that was introduced in Chapter 5.

There are a great many practical applications 41 which demonstrate successful

lighting techniques but no published model exists. This indicates a piecemeal

undisciplined adaptation of the Euclidean based lighting models. Ray tracing

techniques look promising, as does an adaptation of the Torrance- Sparrow

lighting model that was discussed in chapter 5. A good pixel set lighting model
would open the avenue of complex terrain modeling to a much wider audience.

4
Nature's fractal map?
As evidenced by the fractal pictures that have been published.

100
 3. Fractal Music

In [Ref. 7] Voss demonstrates the application of fractal recursive

techniques to — noise and has produced interesting if not pleasing tonal results.

It is safe to surmise that sound is a roughly textured physical phenomenon and

that it may be possible to create or decipher sound using a fractal model. Such a

discovery would aid science in the area of (rapid ) speech recognition.

4. Fractal Computer Graphics Architectures

It is clear from our discussion that new computer architectures need to
be developed to support the pixel set paradigm and the computational aspects of

fractal functions. Such special architectures require parallel processing capabilities

coupled with vast memory resources. A real-time fractal terrain image generator

is one such architectural possibility.

5. A Better Fractal Definition

Fractal geometry is currently attaining a wide audience. Because of that.

it is time that trained mathematicians tackle the problems associated with the
42
imprecise and unworkable current definition of fractal sets . That definition uses

competing definitions of dimension, each of which is somewhat difficult. A new

definition could be based on a fractal set's functional or statistical qualities. Such

a definition scheme must provide tools to further its workability.

42
Sadly, there has been little attention from the mathematical community, although that is

changing. It is with great timidity that ones accepts fractal geometry without such scrutin>

101
B. CONCLUSIONS
Fractal geometry is an old idea that has found a new application with the

advent of computer imaging techniques. Its acceptance, has spawned a great deal

of research and has provided a new tool to observe nature through a different

perspective. We must be careful to insure that our findings are in fact valid. We
also must begin the coalescence of the many techniques that have been developed

in order to control the growth of this concept and to attain true scientific

acceptance. Without this acceptance the theory will be criticized (validly) as an

imprecise and unproven idea . This would be an unfortunate occurrence because

of the potential that fractal geometry possesses.

It is the hope of the author that this work has illuminated the subject of

fractal geometry and that it will aid others in their research. The purpose and
essence of fractal geometry is based on simple concepts. The reader must not be

overawed by the current literature and should retain his perspective with a mild

dose of skepticism. He must not be blinded by skepticism though as the potential

of fractal geometry has not yet been realized. In the final analysis, we expect
that even the skeptical reader will discover the mathematical beauty and

applicative power that fractal geometry possesses.

This of course is the current state of affairs with fractal geometry.

102
 APPENDIX A: FRACTAL COMPUTATION IN R 2

The first routine is the main routine which initializes the data for the Koch curve
generator and initiates the recursive process on each side of the initiator triangle.

The second routine is the recursive subroutine which performs the generator

replacement until the recursive termination event is reached. The termination

event is defined by the precision of the desired output medium.

KOCH.C
/*
This is the main program which controls the initialization of
the koch generator parameters and initiates recursive operations
on each side of the initiator triangle.
*/

/* Global generator and initiator data */

int Generator joints;
/* The number of points in the GENERATOR */
double Gen_angle[lO];
/* The angle formed between init joint 1

and gen joint */

double Gen_ratio[lO];
/* The between init joint 1 to gen joint and
gen joint to init _point2 */
double Tan_theda[lO];
/* The tangent of the angle formed between
init joint 1 and gen joint */
double Curjoint[20][2];
/* Vertices of initiator structure */
int Object joints jmb;
+
/ The number of vertices of the initiating structure */

^include <math.h> /* Standard UNIX include file for math library */

f define x
$ define y 1

103
 2

/ BEGIN MAIN PROGRAM 7

mainQ
I

/* Local variables */
int I;

Initialize global variables "/

f*

/* Initial points ofthe INITIATORS for demo */

Cur point [0][x] - 4.0;

Cur point [0][y] - 3.0 + sqrt(3.0);

( 'ur point [
1 J
[x] 5.0;
Cur point[l][y] - 3.0;

Cur point [2] [x| - 3.0;

Cur point [2] [y] - 3.0;

4
/* Remember to close the side of the triangle /

Cur point[3][x] - 4.0;

Cur point [3][y] - 3.0 t sqrt(3.0);

Object puts iimb 3;

Generator points = 3;

/' Angle (in radians formed between init point 1 and gen point
for demo) '
/

Gen angle[l] 0.0;

Gen angle[2] 0.4712388;

Gen angle[3] 0.0;

I* Ratio of distance between init point and gen point(i) and 1

distance between gen point (i) and init point

( icn rat iol l] 0.5;

Gen ratio[2] -
1.0;

Gen ratiof3 2.0:

it) 4
/<•
r '
'

"
i

r v
 GENERATE.C
/*
This subroutine computes the generator from a given set
of points in R 2
that define a line segment which is the
initiator. The routine is recursive and terminates at a predefined
precision that is input to the subroutine.
7
/* External global generator data; defined in main subroutine */
extern int Generator joints;
/* The number of points in the GENERATOR */
extern double Gen_angle[lO|;
/* The angle formed between init joint 1
and gen joint */
extern double Gen_ratio[lO];
/* The between init joint 1 to gen joint and
gen joint to init joint2 */
extern double Tan_theda[lO];
/* The tangent of the angle formed between
init joint 1 and gen joint */

^include <math.h> /* Standard math include file for UNIX lib */

/* BEGIN RECURSIVE PROCESS */

generate(Xl,Yl,X2,Y2.precision)
/* Parameter variables */
double XI, Yl,X2,Y2, precision;
{

/* Local variables */
long N;
double Parray[3][2];
double G_point[lO][2],DIST;
double Slope _init, Slope perp, Slope gen;
double X jerp,Y jerp,b jerp,b_gen,TEMP:
double ten_thousand, one, zero, minus jne;
int I, J;

/* assign constants */
ten thousand = 10000.0; one = 1.0; zero = 0.0; minus one = -1.0:

106
/* The Koch curve is defined in the infinite but our re<-

willterminate after the distance between points becomes

then the length of a pixel. */

/* Determine distance between point 1 and point 2 "/

TEMP = (X2 - X1)*(X2 - XI - Y2

j
- Yl "
Y2 - Yl
DIST = sqrt( TEMP ;

/* IF DIST less than the precision then terminate this

recursion and begin backtracking */

if (DIST < precision)

{
"
/* Put your Point plotting routine here
printf( "polyline 2"):
c
printf(" xf %l 0.000000". Xl.Yl :

c cw
printf(" 7.f cf 0.000000". X2.Y2 :

return;

/* Put INITIATOR points one and two into the first ar. :

points of the GENERATOR points array as they are alw -

part of the generated structure */

Gjoint'l[l] = XI;
G_point[l][2J= Yl;
G_point[Generator_points +- 2j lj = X2:
G_point[Generator_points -2 2= Y2:

107
/•* Determine the slope of the line formed by the init _pointl
and init joint2. This is the slope of the INITIATOR */

if (X2 != XI)
{
if (Y2 != Yl)
{
Slope Jnit = (Y2 - Yl)/(X2 - Xl);
}
else

{
Slope_init = 0.0;

}
else

{
/* We can't have infinity in a register
so settle with 10k */
Slope_init = ten_thousand;
}

/* For each GENERATOR point (except end points as they are equal
to the INITIATOR end points) find the X,Y values. This is

accomplished by using the data from the global external variables.

The constant data about the ratios and angles between the
INITIATOR and GENERATOR remain the same regardless of the
INITIATORS length or position in EUCLIDIAN space */

for (1=1; I <= Generator_points; I+ + )

{

/* Using the ratios of the generator perpendicular intercept

points on the INITIATOR determine the X,Y values of the
point of intersection of the perpendicular from the
GENERATOR point to the INITIATOR line. */

Xjerp = (XI + Gen_ratio[I]

+
X2)/(1.0 + Gen_ratio[I]);
Yjerp = (Yl + Gen_ratio[I] * Y2)/(1.0 + Gen_ratio[lj);

108
/* If the angle of the INITIATOR point 1 and the GENERATOR
point in question is zero then the GENERATOR point is

coincident with the INITIATOR line and no further

calculations are necessary '/

if ( Gen_angle[I] == zero )

G_point[I+l][l] = X_perp:
G_point[I+l][2] - Yjperp:
}
else

/* There are three STATES possible at this time. STATE 1

where the slope of the initiator line is parallel
to the X or Y axis (which causes havoc with the line
equations). STATE 2 where the slope of the line formed
by. the initiator point 1 and the unknown generator point
is parallel to the X or Y axis. Or STATE 3 where no lines

are parallel to any axis. */

/* Determine the slope of the line through the INITIATOR

point and the unknown
1 GENERATOR point using the
tangent of the Gen_angle in Init.h */

Slope_gen = (Tan_theda[I] + Slope_init)/

(one - Tan_theda[I] * Slope_init);

if ((Slope_gen != zero) &&

(Slope_gen < ten_thousand ))

{
/* Condition one of STATE 3 */

/* Determine Y-intercept for the generator line */

b_gen = Yl - (Slope_gen * Xl);

if ((Slope_init == zero) |

(Slope_init == ten_thousand))
{
/* STATE 1 */

109
 if (Slope_init == ten_thousand )

{
/* STATE 1 condition 1; INITIATOR is parallel
Y axis
to the */

G_point[I+l][2] = Yj>erp;
G_point[I+l][l] = (G_point[I+l][2] - b_gen)/
Slope_gen;
}
else

{
/* STATE 1 condition 2; INITIATOR is parallel
X axis
to the */

G_point[I+l][l] = X j>erp;
G_point[I+l][2] = Slope_gen *

G_point[I+l][l] + b_gen;
}

}
/* END STATE 1 */
else

{
/* STATE 3 */

/* Determine slope of perpendicular line through the

INITIATOR perpendicular intercept. */

Slope_perp = (minus_one)/Slope_init;

/* Determine Y-intercept for perpendicular line */

b_perp = Y_perp - (Slope jerp *

X_perp);

/* Determine the X,Y values of the unknown GENERATOR

point. */

G_point[I+l][l] = (bjperp - b_gen )/

(Slope_j;en - Slope _j>erp);
G_point[I+l][2] = Slope_gen
*

G_point[I+l][l] + b_gen;
}

}
/* END STATE 3 cond. 1 if */

110
 else

{
/* STATE 2 */

Slope_perp = (minus_one)/Slope_init;
b_perp = Y^perp - (Slope j>erp * X_perp):

if (Slope_gen == one)
{
G_point[I+l][l] = XI;
G_point[I+l][2] = Slope^perp * G_point[I+l][l]
+ b_perp;
}
else

{
G_point[I+l][2] = Yl;
G_point[I+l][l] = (G_point[I+l][2| - b_perp)/
Slope_perp;
}

} /* END IF 7
} /* END FOR 7
/* Start recursion on each line formed by the generator from
right to left */

for (J=l; J <= Generator_points + 1; J++)

{
generate(Gj?oint[J] [l],G_point[Jj [2],

G_point[J + l][l],G_point[J+l][2]. precision):

}
/* END generate */
}

111
 )

APPENDIX B: RANDOM NUMBER GENERATORS

The routine below is a C UNIX UCB implementation of the

uniform distribution [0,1] -* standard normal [— oo,oo] transformation. It

generates a 500 entry table of random numbers that observes the period of the

standard normal distribution. Following this routine are statistics that verify the

transformation.

RANDOM. TABLE. GENERATOR.C

I*
This subroutine will build a table in memory that contains 500 random
numbers that observe the period of a standard normal variable
7

^include <math.h> /* Standard UNIX include file for math library */

/* External global variables */

extern double RAND [500];

/* BEGIN MAIN PROGRAM */

rand_table_gen (

{
/* Local Variables */

int I,J;

double UNF1, UNF2;

double range, pi;
int factor;

pi = 3.1415926535;

112
/* Determine the range for the random numbers of UNIX UCB "/

range = 2;
for (J=l; J<=30; J++)

{
range = range *
2;

}
range = range - 1;

/* Set the random number generator seed */

srandom(475836);

/* Create a Table for 500 entries */

for (1=0; I< 500; 1 = 1 + 2)

{
/* Get a uniform random number through the Unix C subroutine */

UNF1 = randomQ;
UNF2 = randomQ;

/* Normalize the uniform random number to the interval [0,1] */

UNF1 = UNF1 / range;

UNF2 = UNF2 / range;

/* Mold the uniform random variable into the approximate normal

distribution */

factor = 1.0;
if (log(UNFl) < 0.0) factor = -1.0;
RAND[I] = sqrt(factor
+
(2.0 * log(UNFl))) *

cos ((2*pi*UNF2));
RAND[I] = RAND[I] * factor;
factor = 1.0:

if (log(UNF2) < 0.0) factor = -1.0;

RAND[I+1] = sqrt(factor * (2.0 * log(UXFl))) *

sin ((2*pi*UNF2));
RAXD[I+1] = RAND[I+1] * factor;

}
return;

113
 VERIFYING STATISTICS
The UNIX UCB operating system's uniform distribution random number
generating function spans the interval defined by its integer range. For a VAX
11/780 implementation this is equivalent to 2 31 - 1 or [0,2147483647].

The random number seed was assigned the value of 475836. The UNIX UCB
random number generator with a fixed seed yields a fixed sequence of numbers

returned from the function, uniformly distributed over the range. This yields a

valuable function if the table needs to be reproduced with the same sequence
after transformation.

The table below shows the results of the uniform distribution sequence after

it was squeezed into the interval [0,1]. These results show that the uniform

distribution has an acceptable distribution over its range. The transformation

31
into [0,1] preserves the distribution from the original range ([0,2 - l]).

Analysis of the normalized uniform random numbers

0.0 -* 0.1 = 52
0.1 -> 0.2 = 47
0.2 - 0.3 = 44
0.3 -> 0.4 = 49
0.4 -+ 0.5 = 49
0.5 -* 0.6 = 47
0.6 -> 0.7 = 51
0.7 -* 0.8 = 57
0.8 - 0.9 = 48
0.9 -> 1.0 = 56

The table below shows the distribution after the

uniform distribution [0,1] -* standard normal [— oo,oc] transformation given the

numbers as described in the above table. This is the data which was used to

build Figure 5.9. The transformation is acceptable for the purpose intended, that

is, to simulate nature's perceived disorder in a fractal function.

114
Analysis of the normal (Gaussian) random numbers

X <= -2.75 = 5
-2.75 < X <= -2.25 = 8
-2.25 < X <= -1.75 = 16
-1.75 < X <= -1.25 = 29
-1.25 < X <= -0.75 = 56
-0.75 < X <= -0.25 = 88
-0.25 < X <= 0.25 = 115
0.25 < X <= 0.75 = 87
0.75 < X <= 1.25 = 59
1.25 < X <= 1.75 = 21
1.75 < X <= 2.25 = 12
2.25 < X <= 2.75 = 4
2.75 < X =0

115
 C
APPENDIX C: THE TRIANGULAR MOUNTAIN

The first routine is the main routine which initializes the generator data for

the initiating triangle and initiates the recursive process. The second routine is

the recursive subroutine which performs the generator replacement until the

recursive termination event is reached, which is defined by the precision

parameter.

MOUNTAIN.
/*
This is the main program that controls the initialization of the triangle
initiating structure and initiates the recursion on that triangle. The
recursion will proceed until the recursive termination event (defined
by the precision global parameter)
*/

^include <math.h> /* Standard UNIX include file for math library */

/* Global Defines */

# define x
f define y 1

^define z 2

/* Global Tables */

double RAND[500];
double Precision;
double Scale;

116
/* BEGIN MAIN PROGRAM */

mainQ

{
/* Local Variables */

int I,J,K;
double P1[3],P2[3],P3[3];
int Seed;

/* Create the initiating triangle structure */

Pl[x] = 4.5;
Pl[y] - 3.25;
Pljz] = 0.0;

P2[x] = 7.0;
P2[y] = 3.25;
P2[z] - 0.0;

P3[x] = 5.75;
P3[y] = 3.25 + sqrt(((2.5 * 2.5) - (1.25
+
1.25)));
P3[z] = 0.0;

/* Build the random number table (appendix B) */

rand_table_gen();

/* Fractalize until desired precision */

Seed = 0; /* Entry seed to the random number table */

Precision = 0.3; /* Recursive termination distance */
Scale = 0.2; /* Scaling factor for vertical Y displacement
in the mountain_generate subroutine */

mountain_generate(Pl,P2.P3.Seed);

/* END MAIN */

117
 ; C

GENERATE MOUNTAIN.
/*
This is the subroutine that computes the four generated triangles from
an initiating triangle. The routine is recursive and terminates at a
predefined precision defined in the global parameter Precision
r
^include <math.h> /* Standard math include file for UNIX lib */

f define x
# define y 1

# define z 2

/* Global Structures */

extern double RAND[500];

extern double Precision;
extern double Scale;

/* BEGIN RECURSIVE PROCESS */

mount ain_generate(Pl,P2,P3, Seed)

/* Parameter variables */

double P1[3],P2[3],P3[4];
int Seed;

/* Local variables */

int I,J;

double Pgenl [3] ,Pgen2 [3] ,Pgen3 [3]

double Pmidl[3],Pmid2[3].Pmid3[3];
double TEMP,DIST,TWO;

TWO = 2.0;

118
 ; /

/* Determine distance between point 1 and point 2 *

TEMP = (P2[xj - Pl[x])

+
(P2[x] - Pl[x|) +
(P2[y|-Pl[y])*(P2[y]-Pl[y]) +
(P2[z] - Pl[z])*(P2[z] - Pl[z]);
DIST = sqrt( TEMP );

/* IF DIST less than one then terminate this recursion and

begin backtracking */

if (DIST < Precision)

{
/* Put your polygon plotting routine here */
printf "poly gon3 "
(
)

MO/
printf( of %f %f\Pl[x],Pl[y],Pl[z]):
printf("%f %f %f \P2[x],P2[y],P2[zj):
printf("%f %f %P',P3[x],P3[y],P3[z]);
return;

/* Manage the Seed number for a 500 entry table */

if (Seed > 496) Seed = 0;

/* Find the midpoints of each triangle leg */

for (1=0; I<=2; I+ + ) /* thru 2 => x.y,z */

{
Pmidl[I] = (P1[I] + P2[I]) / TWO;
}

for (1=0; I<=2; I+ + )

/* thru 2 => x,y,z */

{
Pmid2[I] = (P2[I] + P3[I|) / TWO;
}

for (1=0; I< = 2: I+ + )

/* thru 2 => x.v.z */

{
Pmid3[I] = (P3[I] + P1[I]) / TWO;
}

119
/* Adjust the Y coordinate => normal from Z-X plane */

Pmidl[y] = (Scale * RAND[Seed]) + Pmidl[y];

Pmid2[y] = (Scale * RAND[Seed+l]) + Pmid2[y];
Pmid3[y] = (Scale * RAND[Seed+2j) + Pmid3[yj;
Seed = Seed 4- 1;

/* Recurse on the triangles according to the reverse order rule

for the interior triangle to preserve seed order */

mountain_generate(Pmidl,P2, Pmid2,Seed);
mountain_generate(Pmid3,Pmid2,P3. Seed);
mountain_generate(Pl, Pmidl,Pmid3,Seed);
mountain_generate(Pmid2,Pmid3,Pmidl,Seed);
/* END generate */
}

120
 C

APPENDIX D: THE RECTANGULAR MOUNTAIN

The first routine is the main routine which initializes the generator data for

the initiating rectangular shape and initiates the recursive process. The second
routine is the recursive subroutine which performs the generator replacement

until the recursive termination event is reached, which is defined by the precision

parameter.

RECTANGULAR-MOUNTAIN.
/*
This is the main program that controls the initialization of the rectangular
initiating structure and initiates the recursion on that rectangle. The
recursion will proceed until the recursive termination event (defined
by the precision global parameter)
7
^include <math.h> /* Standard UNIX include file for math library */

/* Global Defines */

fdefine x
# define y 1

# define z 2

/* Global Tables */

double RAND [500];

double Precision;
double Scale;

121
/* BEGIN MAIN PROGRAM */

main()
{
/* Local Variables */
int I,J,K;
double P1[3],P2[3],P3[3],P4[3];
int Seed;

/* Create the four sided polygon initiating structure */

Pl[x] = 2.0;
Pl[y] = 5.0;
Pljz] = 1.0;

P2[x] = 6.0;
P2[y] = 5.0;
P2[z] = 1.0;

P3[x] = 7.5;
P3[y] = 6.5;
P3[z] = 3.0;

P4[x] = 3.5;
P4[y] = 6.5;
P4[z - 3.0;

/* Build the random number table (appendix B) */

rand_table_gen();

/* Fractalize until desired precision */

Seed = 100;
Precision = 0.5;
Scale = 0.07;

mountain_generate(Pl,P2,P3,P4.Seed);

/* END MAIN */

122
 C

RECTANGULAR. GENERATE. MOUNTAIN.

/*
This is the subroutine that computes the four generated rectangles from
an initiating rectangle. The routine is recursive and terminates at a
predefined precision defined in the global parameter Precision
l*

^include <math.h> /* Standard math include file for UNIX lib */

§ define x
$ define y 1
# define z 2

/* Global Structures */

extern double RAND [500];

extern double Precision:
extern double Scale;

/* BEGIN RECURSIVE PROCESS */

mountain_generate(Pl,P2,P3,P4,Seed)

/* Parameter variables */
double P1[3],P2[3].P3[3].P4[3);
int Seed:
{
/* Local variables */
int I,J;

double Pmidl[3].Pmid2[3],Pmid3[3].Pmid4[3],Center[3];
double TEMP,DIST.T\VO.FOUR:

TWO = 2.0: FOUR - 4.0:

/* Determine distance between point 1 and point 2 */

TEMP = (P2[x] - Pl[x)nP2[x] - Pl[x]j +
(P2[y] - Pl[y]r(P2[y] - Pl(y]) +
(P2[z] - Pl[z])*(P2[z] - Pl[z]);
DIST = sqrt( TEMP );

123
 ;

/* If DIST less than one then terminate this recursion and

begin backtracking */
if (DIST < Precision)

{
/* Put your Polygon output routine here */
printf ( "polygon4 " )
printf("%f %f %f",Pl[x],Pl[y],Pl[z]);
printf("%f %f %f\P2[x],P2[y],P2[zj):
printf("%f %f %f\P3[x],P3[y],P3[zj):
printf("%f %f %f*,P4[x],P4[y],P4[z]):
return;

/* Manage the Seed number for a 500 entry table */

if (Seed > 496) Seed = 0;

/* Find the midpoints of each rectangle leg */

for (1=0; I<=2; I+ + /* thru 2 => x,y,z */
)

{
Pmidl[I] = (P1[I] + P2[I]) / TWO;
}

for (1=0; I< = 2; I+ + )

/* thru 2 => x,y,z */

{
Pmid2[I] = (P2[I] + P3[I]) / TWO;
}

for (1=0; I<=2; I+ + ) /* thru 2 => x,y,z */

{
Pmid3[I] = (P3[I] + P4[I]) / TWO;
}

for (1=0; I<=2; I+ + ) /* thru 2 => x.y,z */

{
Pmid4[I] = (P1[I] + P4[I]) / TWO:
}

124
/* The four sided polygon is non-planar so average the xyz-displacement
for a best approach */
fit

for (1=0; I<=2; I+ + ) /* thru 2 => x.v.z '/

{
Center[I] = (P3[I] + PI [I] + P4[I] + P2[I]) / FOUR;
}

/* Adjust the Y coordinate => normal from Z-X plane */

Pmidl[y] = (Scale * RANDjSeed]) + Pmidl[y);
Pmid2[y] = (Scale * RAND[Seed+l]) + Pmid2[y];
Pmid3[y] = (Scale * RAND[Seed+2]) + Pmid3(yj;
Pmid4[y] = (Scale * RAND[Seed+3]) + Pmid4[yj;
Center[y] = (Scale * RAND[Seed+4]) + Center[y];
Seed = Seed + 4;

/* Recurse on the rectangles according to the reverse order rule

for the interior rectangles to preserve seed order */

mountain_generate(Pl, Pmidl. Center, Pmid4, Seed)

mountain_generate(Pmid2, Center, Pmidl, P2, Seed)
mountain_generate(Pmid2, Center, Pmid3, P3, Seed)
mountain_generate(P4, Pmid3, Center, Pmid4, Seed)

/* END generate */
}

125
 APPENDIX E: GEOMETRIC SUPPORT
Many fractal applications and computer graphics models use the normal to a

plane as a computational reference point. For this reason, this appendix is

devoted to two tools for determining the plane equation of a polygon and the

equation of the normal to the computed plane.

Determinant Approach to the Planar Equation

One of the most common forms of a planar equation is the general form. This

form uniquely describes a plane through four coefficients A,B,C and D:

Ax + By + Cz = D

With three points on a plane, you can determine the planar equation by
computing the coefficients. This approach utilizes the determinant form of the

planar equation. Given the points Pj = (x l ,y 1 ,z x ), P2 = ( j 2 -?/ 2->-s 2) an<^

P3 = (£3,2/3,23) such that Pj ^ P2 # P3, these points determine a unique plane

in space through the determinant equation:

x - x1 y - y l
z - zj

x2 ~ x l Y2 ~ Yl z2 ~ zl =
x2 - x i
y~2 - yi z2 - zi

To simplify the equation, we replace the constant differences by the expressions:

Cl x - x 2 —
C2 X = x 3 -
Cl y = x 2 -
C2 y = x 3 -
CI. = x 9 -

C2. = Xo -

126
Evaluating the determinant using the diagonal approach yields:

[(x - x 1 )Cl C2, - (x - x,)Cl z C2 +

y y ]

[(y - yi )Cl z C2 x - (y - yJCl»C2j +

[(z - Zl )Cl x C2 y - (z - Bl )Cl y C2j =

Solving the equations for x,y and z in terms of the constant expressions:

A = Cl C2 z - Cl z C2
y y

B = Cl z C2 x - C1 X C2 Z
C = Cl x C2 y - Cl y C2 x
D = -[A Xl + B yi + CzJ

The Normal to the Plane

Once the parameters A,B and C have been determined, the solution of the

linear equation for any normal to the plane is straightforward. Using the plane

parameters in the parametric equation for the normal line to the plane and using

any known point on the plane (x kwn ,y kwn .z kwn )

(the midpoint of the fractal

triangle for example) determines a normal line as:

X = * kwn + c A
Y = y kwn + c B
Z = z kwn + c C

where c is a parameter such that c is an element of R. By varying the parameter

c we can solve for unique points on the normal line to the plane.

127
 LIST OF REFERENCES

1. Benoit B. Mandelbrot, The Fractal Geometry of Nature, W. H.

Freeman and Company, 1983.

2. Donald F. Stanet and David F. McAllister Discrete Mathematics

in Computer Science, Prentice-Hall, Inc., 1977.

3. Felix Hausdorff, Set Theory, Chelsea Publishing Company, 1957.

4. Alan Norton, Generation and Display of Geometric Fractals in 3-D.

Computer Graphics, vol 18, no. 13, July 1984.

5. J. Perrin, La Discontinuite s de la Matiere, Revue du Mois, vol 1, 1906;

quoted by Mandelbrot in Ref. l:pp 7-9.

6. Witold Hurewicz and Henry Wallman, Dimension Theory, Princeton

University Press, 1941.

7. Special Interest Group on Computer Graphics of the Association for

Computing Machinery (SIGGRAPH), Fractals: Basic Concepts,
Computation and Rendering, Course on Fractals July 23, 1985.

8. James D. Foley and Andries Van Dam, Fundamentals of Interactive

Computer Graphics, Addison- Wesley, 1982.

9. Jagdish K. Patel and Campbell B. Read, Handbook of the Normal

Distribution, Marcel Dekker, Inc, 1982.

128
 INITIAL DISTRIBUTION LIST

No. Copies

1. Defense Technical Information Center 2

Cameron Station
Alexandria, Virginia 22304-6145

2. Superintendent 2
Attn: Library (Code 0142)
Naval Postgraduate School
Monterey, California 93943

3. Chairman (Code 52) 1

Department of Computer Science

Naval Postgraduate School
Monterey, California 93943

4. Computer Technology Curricular Officer (Code 37) 1

Naval Postgraduate School

Monterey, California 93943

*
5. Michael Zyda (Code 52Zk)
J. 2
Department of Computer Science
Naval Postgraduate School
Monterey, California 93943

6. Daniel Davis (Code 52Dv) 1

Department of Computer Science

Naval Postgraduate School
Monterey, California 93943

7. Capt. Michael E. Gaddis 10

DCA/JDSSC/C720
1860 Wiehle Ave.
Reston Virginia 22090

8. Lt. John M. Yurchak 1

P.O. Box 268 RD#4

Muncy, Pennsylvania 17756

9. Capt. H. McArthur 1

11243 Edgemoor Ct.

Woodbridge Virginia 22192

129
10. Richard W. Hamming (Code 52Hg)
Department of Computer Science
Naval Postgraduate School
Monterey, California 93943

11. Bruce J. MacLennan (Code 52M1)

Department of Computer Science
Naval Postgraduate School
Monterey, California 93943

12. Vincent Lum (Code 52V1)

Department of Computer Science
Naval Postgraduate School
Monterey, California 93943

13. Lt(N) Peter J. Lenk

DMCS 3-2-6
101 Colonel By Dr.
Ottawa, Ontario, Canada

14. Leut. Steven J. Firth RAN

15 Carrington Rd.
Hornsby, New South Wales, Australia

15. Robert B. McGhee (Code 52Mz)

Department of Computer Science
Naval Postgraduate School
Monterey, California 93943

16. Larry J. Williamson (Code 53Wj)

Department of Mathematics
Naval Postgraduate School
Monterey, California 93943

130
 Thesis
G1223
216443
Th'c.l Gaddis
iGl
of nature- e
7 S °metr y
I
m tica i

Vr
the -
f basis ?n
plic ap ~
ation to
Com Puter
graphics.

lft 02g

216443
Gaddis
The fractal geometry
of nature; its mathe-
matical basis and ap-
plication to computer
graphics.
thesG1223
The fractal geometry of nature; its mat

3 2768 000 65128 5

DUDLEY KNOX LIBRARY