Aristid Lindenmayer 1925–1989

Przemyslaw PrusinkiewiczAristid Lindenmayer

The Algorithmic Beautyof Plants

WithJames S. HananF. David FracchiaDeborah FowlerMartin J. M. de BoerLynn Mercer

With 150 Illustrations, 48 in Color

This edition of The Alogirthmic Beauty of Plants is the electronic versionof the book that was published by Springer-Verlag, New York, in 1990 andreprinted in 1996. The electronic version has been produced using the orig-inal LATEX files and digital illustrations.

c©2004 Przemyslaw PrusinkiewiczAll rights reserved.

Front cover design: The roses in the foreground (Roses by D. R. Fowler,J. Hanan and P. Prusinkiewicz [1990]) were modeled using L-systems. Dis-tributed ray-tracing with one extended light source was used to simulatedepth of field. The roses were placed on a background image (photgraphyby G. Rossbach), which was scanned digitally and post-processed.

Preface

The beauty of plants has attracted the attention of mathematicians for Mathematicsand beautycenturies. Conspicuous geometric features such as the bilateral sym-

metry of leaves, the rotational symmetry of flowers, and the helicalarrangements of scales in pine cones have been studied most exten-sively. This focus is reflected in a quotation from Weyl [159, page 3],“Beauty is bound up with symmetry.”

This book explores two other factors that organize plant structuresand therefore contribute to their beauty. The first is the elegance andrelative simplicity of developmental algorithms, that is, the rules whichdescribe plant development in time. The second is self-similarity, char-acterized by Mandelbrot [95, page 34] as follows:

When each piece of a shape is geometrically similar to thewhole, both the shape and the cascade that generate it arecalled self-similar.

This corresponds with the biological phenomenon described by Herman,Lindenmayer and Rozenberg [61]:

In many growth processes of living organisms, especially ofplants, regularly repeated appearances of certain multicel-lular structures are readily noticeable.... In the case of acompound leaf, for instance, some of the lobes (or leaflets),which are parts of a leaf at an advanced stage, have thesame shape as the whole leaf has at an earlier stage.

Thus, self-similarity in plants is a result of developmental processes. Growth andformBy emphasizing the relationship between growth and form, this book

follows a long tradition in biology. D’Arcy Thompson [143] traces itsorigins to the late seventeenth century, and comments:

Organic form itself is found, mathematically speaking, to bea function of time.... We might call the form of an organisman event in space-time, and not merely a configuration inspace.

This concept is echoed by Halle, Oldeman and Tomlinson [58]:

The idea of the form implicitly contains also the history ofsuch a form.

vi Preface

The developmental processes are captured using the formalism ofL-systems. They were introduced in 1968 by Lindenmayer [82] as aModeling of

plants theoretical framework for studying the development of simple multicel-lular organisms, and subsequently applied to investigate higher plantsand plant organs. After the incorporation of geometric features, plantmodels expressed using L-systems became detailed enough to allow theuse of computer graphics for realistic visualization of plant structuresand developmental processes.

The emphasis on graphics has several motivations. A visual compar-ison of models with real structures is an important component of modelvalidation. The display of parameters and processes not observable di-rectly in living organisms may assist in the analysis of their physiology,and thus present a valuable tool for educational purposes. From anaesthetic perspective, plants present a wealth of magnificent objectsfor image synthesis. The quest for photorealism challenges modelingand rendering algorithms, while a departure from realism may offer afresh view of known structures.

The application of computer graphics to biological structures is onlyone of many factors that contribute to the interdisciplinary characterof this book. For example, the notion of L-systems is a part of formallanguage theory, rooted in the theory of algorithms. The application ofL-systems to plant description has been studied by biologists, and in-volves various methods of general mathematics. Self-similarity relatesplant structures to the geometry of fractals. Computer-aided visual-ization of these structures, and the processes that create them, joinsscience with art.

The study of an area that combines so many disciplines is very stim-About the bookulating. Some results may be of special interest to students of biologyor computer graphics, but a much wider circle of readers, generally in-terested in science, may find mathematical plant models inspiring, andthe open problems worth further thought. Consequently, all basic con-cepts are presented in a self-contained manner, assuming only generalknowledge of mathematics at the junior college level.

This book focuses on original research results obtained by the au-thors in the scope of the cooperation between the Theoretical BiologyGroup, directed by Aristid Lindenmayer at the University of Utrecht,and the Computer Graphics Group, working under the supervision ofPrzemyslaw Prusinkiewicz at the University of Regina. Technically, thebook evolved from the SIGGRAPH ’88 and ’89 course notes Linden-mayer systems, fractals, and plants, published by Springer-Verlag inthe series Lecture Notes in Biomathematics [112]. The present volumehas been extended with edited versions of recent journal and conferencepapers (see Sources), as well as previously unpublished results.

Aristid Lindenmayer is the author of the notion of L-systems whichforms the main thread of the book. He also played an essential rolein the reported research by suggesting topics for study, guiding theconstruction of specific plant models, monitoring their correctness and

Preface vii

participating in many discussions of biological and mathematical prob-lems. Seriously ill, Professor Lindenmayer co-authored and edited sev-eral chapters, but was not able to participate in the completion of thiswork. If any inaccuracies or mistakes remain, he could not preventthem. Still, in spite of unavoidable shortcomings, we hope that thisbook will convey his and our excitement of applying mathematics toexplore the beauty of plants.

viii Preface

Acknowledgements

While preparing this book, we received extraordinary support and helpfrom many people, and we are deeply thankful to all of them. Firstof all, we would like to thank those who were directly involved in theunderlying research and software development. Craig Kolb wrote theray tracer rayshade used to render many of the images included inthe book. Allan Snider developed several software tools, including apreviewer for rayshade, and provided valuable expertise in ray-tracing.Daryl Hepting developed software for rendering sets defined by iteratedfunction systems and provided diagrams for Chapter 8. Norma Fullermodeled several man-made objects incorporated into the images.

We would like to thank Zdzis�law Pawlak and Grzegorz Rozenbergwho initiated the contact between the Theoretical Biology Group atthe University of Utrecht and the Computer Graphics Group at theUniversity of Regina. Benoit Mandelbrot and Heinz-Otto Peitgen madeit possible to conduct parts of the reported research at Yale Universityand the University of Bremen.

We are also grateful to all those who shared their knowledge withus and made suggestions reflected in this book. Discussions and corre-spondence with Jules Bloomenthal, Mark de Does, Pauline Hogeweg,Jacqueline and Hermann Luck, Gavin Miller, Laurie Reuter, DietmarSaupe and Alvy Ray Smith were particularly fruitful.

Research reported in this book was funded by grants from the Nat-ural Sciences and Engineering Research Council of Canada, as well asan equipment donation and a research grant from Apple Computer,Inc. We are particularly grateful to Mark Cutter for making the sup-port from Apple possible. The Graphics Laboratory at the Universityof Regina also enjoys continued support from the university. The influ-ence of Lawrence Symes and R. Brien Maguire is deeply appreciated.In addition, the University of Regina and the University of Utrechtcontributed towards travel expenses.

We would like to thank Springer-Verlag and in particular, GerhardRossbach and Nina LaVoy from the Springer West Coast Office, for theexpedient publishing of this book.

Finally, we would like to thank our families and friends for theirlove, support and patience while we worked on this book.

Przemyslaw PrusinkiewiczJames HananF. David FracchiaDeborah R. FowlerMartin J. M. de BoerLynn MercerRegina, CanadaMay 1990

Sources ix

Sources

In the preparation of this book, edited parts of the following publica-tions were used:

• P. Prusinkiewicz and J. Hanan. Lindenmayer systems, fractals,and plants, volume 79 of Lecture Notes in Biomathematics.Springer-Verlag, Berlin, 1989. [Chapters 1, 3, 5]

• P. Prusinkiewicz, A. Lindenmayer and J. Hanan. Developmen-tal models of herbaceous plants for computer imagery purposes.Proceedings of SIGGRAPH ’88 (Atlanta, Georgia, August 1-5,1988), in Computer Graphics, 22(4):141–150, 1988. Used with thepermission of the Association for Computing Machinery. [Chap-ters 1, 3, 5]

• P. Prusinkiewicz and J. Hanan. Visualization of botanical struc-tures and processes using parametric L-systems. In D. Thal-mann, editor, Scientific visualization and graphics simulation,pages 183–201. J. Wiley & Sons, 1990. [Chapters 1–5]

• D. Fowler, J. Hanan and P. Prusinkiewicz. Modelling spiral phyl-lotaxis. computers & graphics, 13(3):291–296, 1989. [Chapter 4]

• J. S. Hanan. PLANTWORKS: A software system for realisticplant modelling. Master’s thesis, University of Regina, 1988.[Chapter 5]

• F. D. Fracchia, P. Prusinkiewicz and M. J. M. de Boer. Animationof the development of multicellular structures. In N. Magnenat-Thalmann and D. Thalmann, editors, Computer Animation ’90,pages 3–18. Springer-Verlag, Tokyo, 1990. [Chapter 7]

• F. D. Fracchia, P. Prusinkiewicz, and M. J. M. de Boer. Visualiza-tion of the development of multicellular structures. In Proceedingsof Graphics Interface ’90, pages 267–277, 1990. [Chapter 7]

• L. Mercer, P. Prusinkiewicz, and J. Hanan. The concept and de-sign of a virtual laboratory. In Proceedings of Graphics Interface’90, pages 149–155, 1990. [Appendix 8.2]

Contents

1 Graphical modeling using L-systems 11.1 Rewriting systems . . . . . . . . . . . . . . . . . . . . . . 11.2 DOL-systems . . . . . . . . . . . . . . . . . . . . . . . . 31.3 Turtle interpretation of strings . . . . . . . . . . . . . . . 61.4 Synthesis of DOL-systems . . . . . . . . . . . . . . . . . 11

1.4.1 Edge rewriting . . . . . . . . . . . . . . . . . . . 111.4.2 Node rewriting . . . . . . . . . . . . . . . . . . . 131.4.3 Relationship between edge and

node rewriting . . . . . . . . . . . . . . . . . . . . 181.5 Modeling in three dimensions . . . . . . . . . . . . . . . 181.6 Branching structures . . . . . . . . . . . . . . . . . . . . 21

1.6.1 Axial trees . . . . . . . . . . . . . . . . . . . . . . 211.6.2 Tree OL-systems . . . . . . . . . . . . . . . . . . 231.6.3 Bracketed OL-systems . . . . . . . . . . . . . . . 24

1.7 Stochastic L-systems . . . . . . . . . . . . . . . . . . . . 281.8 Context-sensitive L-systems . . . . . . . . . . . . . . . . 301.9 Growth functions . . . . . . . . . . . . . . . . . . . . . . 361.10 Parametric L-systems . . . . . . . . . . . . . . . . . . . . 40

1.10.1 Parametric OL-systems . . . . . . . . . . . . . . . 411.10.2 Parametric 2L-systems . . . . . . . . . . . . . . . 431.10.3 Turtle interpretation of parametric words . . . . . 46

2 Modeling of trees 51

3 Developmental models of herbaceous plants 633.1 Levels of model specification . . . . . . . . . . . . . . . . 64

3.1.1 Partial L-systems . . . . . . . . . . . . . . . . . . 643.1.2 Control mechanisms in plants . . . . . . . . . . . 653.1.3 Complete models . . . . . . . . . . . . . . . . . . 68

3.2 Branching patterns . . . . . . . . . . . . . . . . . . . . . 703.3 Models of inflorescences . . . . . . . . . . . . . . . . . . 71

3.3.1 Monopodial inflorescences . . . . . . . . . . . . . 713.3.2 Sympodial inflorescences . . . . . . . . . . . . . . 823.3.3 Polypodial inflorescences . . . . . . . . . . . . . . 863.3.4 Modified racemes . . . . . . . . . . . . . . . . . . 93

xii Contents

4 Phyllotaxis 994.1 The planar model . . . . . . . . . . . . . . . . . . . . . . 1004.2 The cylindrical model . . . . . . . . . . . . . . . . . . . . 109

5 Models of plant organs 1195.1 Predefined surfaces . . . . . . . . . . . . . . . . . . . . . 1195.2 Developmental surface models . . . . . . . . . . . . . . . 1205.3 Models of compound leaves . . . . . . . . . . . . . . . . 128

6 Animation of plant development 1336.1 Timed DOL-systems . . . . . . . . . . . . . . . . . . . . 1356.2 Selection of growth functions . . . . . . . . . . . . . . . . 139

6.2.1 Development of nonbranching filaments . . . . . . 1406.2.2 Development of branching structures . . . . . . . 142

7 Modeling of cellular layers 1457.1 Map L-systems . . . . . . . . . . . . . . . . . . . . . . . 1457.2 Graphical interpretation of maps . . . . . . . . . . . . . 1507.3 Microsorium linguaeforme . . . . . . . . . . . . . . . . . 1537.4 Dryopteris thelypteris . . . . . . . . . . . . . . . . . . . . 1627.5 Modeling spherical cell layers . . . . . . . . . . . . . . . 1667.6 Modeling 3D cellular structures . . . . . . . . . . . . . . 168

8 Fractal properties of plants 1758.1 Symmetry and self-similarity . . . . . . . . . . . . . . . . 1778.2 Plant models and iterated function systems . . . . . . . . 178

Epilogue 191

Appendix A Software environment for plant modeling 193A.1 A virtual laboratory in botany . . . . . . . . . . . . . . . 193A.2 List of laboratory programs . . . . . . . . . . . . . . . . 198

Appendix B About the figures 201

Appendix C Turtle interpretation of symbols 209

Bibliography 211

Index 225

Chapter 1

Graphical modeling usingL-systems

Lindenmayer systems — or L-systems for short — were conceived asa mathematical theory of plant development [82]. Originally, they didnot include enough detail to allow for comprehensive modeling of higherplants. The emphasis was on plant topology, that is, the neighborhoodrelations between cells or larger plant modules. Their geometric aspectswere beyond the scope of the theory. Subsequently, several geometricinterpretations of L-systems were proposed with a view to turning theminto a versatile tool for plant modeling. Throughout this book, aninterpretation based on turtle geometry is used [109]. Basic notionsrelated to L-system theory and their turtle interpretation are presentedbelow.

1.1 Rewriting systems

The central concept of L-systems is that of rewriting. In general, rewrit-ing is a technique for defining complex objects by successively replacingparts of a simple initial object using a set of rewriting rules or produc-tions . The classic example of a graphical object defined in terms ofrewriting rules is the snowflake curve (Figure 1.1), proposed in 1905 by Koch

constructionvon Koch [155]. Mandelbrot [95, page 39] restates this construction asfollows:

One begins with two shapes , an initiator and a generator.The latter is an oriented broken line made up of N equalsides of length r. Thus each stage of the construction beginswith a broken line and consists in replacing each straightinterval with a copy of the generator, reduced and displacedso as to have the same end points as those of the intervalbeing replaced.

2 Chapter 1. Graphical modeling using L-systems

initiator

generator

Figure 1.1: Construction of the snowflake curve

While the Koch construction recursively replaces open polygons, rewrit-ing systems that operate on other objects have also been investigated.For example, Wolfram [160, 161] studied patterns generated by rewrit-ing elements of rectangular arrays. A similar array-rewriting mecha-nism is the cornerstone of Conway’s popular game of life [49, 50]. Animportant body of research has been devoted to various graph-rewritingsystems [14, 33, 34].

The most extensively studied and the best understood rewriting sys-Grammarstems operate on character strings. The first formal definition of such asystem was given at the beginning of this century by Thue [128], buta wide interest in string rewriting was spawned in the late 1950s byChomsky’s work on formal grammars [13]. He applied the concept ofrewriting to describe the syntactic features of natural languages. Afew years later Backus and Naur introduced a rewriting-based notationin order to provide a formal definition of the programming languageALGOL-60 [5, 103]. The equivalence of the Backus-Naur form (BNF)and the context-free class of Chomsky grammars was soon recognized[52], and a period of fascination with syntax, grammars and their appli-cation to computer science began. At the center of attention were setsof strings — called formal languages — and the methods for generating,recognizing and transforming them.

In 1968 a biologist, Aristid Lindenmayer, introduced a new type ofL-systemsstring-rewriting mechanism, subsequently termed L-systems [82]. Theessential difference between Chomsky grammars and L-systems lies in

1.2. DOL-systems 3

Figure 1.2: Relations between Chomsky classes of languages and languageclasses generated by L-systems. The symbols OL and IL denote languageclasses generated by context-free and context-sensitive L-systems, respec-tively.

the method of applying productions. In Chomsky grammars produc-tions are applied sequentially, whereas in L-systems they are appliedin parallel and simultaneously replace all letters in a given word. Thisdifference reflects the biological motivation of L-systems. Productionsare intended to capture cell divisions in multicellular organisms, wheremany divisions may occur at the same time. Parallel production ap-plication has an essential impact on the formal properties of rewritingsystems. For example, there are languages which can be generatedby context-free L-systems (called OL-systems) but not by context-freeChomsky grammars [62, 128] (Figure 1.2).

1.2 DOL-systems

This section presents the simplest class of L-systems, those which aredeterministic and context-free, called DOL-systems. The discussionstarts with an example that introduces the main idea in intuitive terms.

Consider strings (words) built of two letters a and b, which may Exampleoccur many times in a string. Each letter is associated with a rewritingrule. The rule a → ab means that the letter a is to be replaced bythe string ab, and the rule b → a means that the letter b is to bereplaced by a. The rewriting process starts from a distinguished stringcalled the axiom. Assume that it consists of a single letter b. In thefirst derivation step (the first step of rewriting) the axiom b is replaced

4 Chapter 1. Graphical modeling using L-systems

Figure 1.3: Example of a derivation in a DOL-system

by a using production b → a. In the second step a is replaced by abusing production a → ab. The word ab consists of two letters, both ofwhich are simultaneously replaced in the next derivation step. Thus, ais replaced by ab, b is replaced by a, and the string aba results. In asimilar way, the string aba yields abaab which in turn yields abaababa,then abaababaabaab, and so on (Figure 1.3).

Formal definitions describing DOL-systems and their operation aregiven below. For more details see [62, 127].

Let V denote an alphabet, V ∗ the set of all words over V , andL-systemV + the set of all nonempty words over V . A string OL-system is anordered triplet G = 〈V, ω, P 〉 where V is the alphabet of the system,ω ∈ V + is a nonempty word called the axiom and P ⊂ V × V ∗ is afinite set of productions. A production (a, χ) ∈ P is written as a →χ. The letter a and the word χ are called the predecessor and thesuccessor of this production, respectively. It is assumed that for anyletter a ∈ V , there is at least one word χ ∈ V ∗ such that a → χ. Ifno production is explicitly specified for a given predecessor a ∈ V , theidentity production a → a is assumed to belong to the set of productionsP . An OL-system is deterministic (noted DOL-system) if and only iffor each a ∈ V there is exactly one χ ∈ V ∗ such that a → χ.

Let µ = a1 . . . am be an arbitrary word over V . The word ν =Derivationχ1 . . . χm ∈ V ∗ is directly derived from (or generated by) µ, noted µ ⇒ν, if and only if ai → χi for all i = 1, . . . ,m. A word ν is generated byG in a derivation of length n if there exists a developmental sequence ofwords µ0, µ1, . . . , µn such that µ0 = ω, µn = ν and µ0 ⇒ µ1 ⇒ . . . ⇒µn.

1.2. DOL-systems 5

Figure 1.4: Development of a filament (Anabaena catenula) simulated usinga DOL-system

The following example provides another illustration of the operation of AnabaenaDOL-systems. The formalism is used to simulate the development of afragment of a multicellular filament such as that found in the blue-greenbacteria Anabaena catenula and various algae [25, 84, 99]. The symbolsa and b represent cytological states of the cells (their size and readinessto divide). The subscripts l and r indicate cell polarity, specifying thepositions in which daughter cells of type a and b will be produced. Thedevelopment is described by the following L-system:

ω : ar

p1 : ar → albr

p2 : al → blar

p3 : br → ar

p4 : bl → al

(1.1)

Starting from a single cell ar (the axiom), the following sequence ofwords is generated:

ar

albr

blarar

alalbralbr

blarblararblarar

· · ·

6 Chapter 1. Graphical modeling using L-systems

Under a microscope, the filaments appear as a sequence of cylin-ders of various lengths, with a-type cells longer than b-type cells. Thecorresponding schematic image of filament development is shown inFigure 1.4. Note that due to the discrete nature of L-systems, the con-tinuous growth of cells between subdivisions is not captured by thismodel.

1.3 Turtle interpretation of strings

The geometric interpretation of strings applied to generate schematicimages of Anabaena catenula is a very simple one. Letters of theL-system alphabet are represented graphically as shorter or longer rect-angles with rounded corners. The generated structures are one-dimen-sional chains of rectangles, reflecting the sequence of symbols in thecorresponding strings.

In order to model higher plants, a more sophisticated graphical in-Previousmethods terpretation of L-systems is needed. The first results in this direction

were published in 1974 by Frijters and Lindenmayer [46], and Hogewegand Hesper [64]. In both cases, L-systems were used primarily to de-termine the branching topology of the modeled plants. The geometricaspects, such as the lengths of line segments and the angle values, wereadded in a post-processing phase. The results of Hogeweg and Hesperwere subsequently extended by Smith [136, 137], who demonstrated thepotential of L-systems for realistic image synthesis.

Szilard and Quinton [141] proposed a different approach to L-systeminterpretation in 1979. They concentrated on image representationswith rigorously defined geometry, such as chain coding [43], and showedthat strikingly simple DOL-systems could generate the intriguing, con-voluted curves known today as fractals [95]. These results were sub-sequently extended in several directions. Siromoney and Subrama-nian [135] specified L-systems which generate classic space-filling curves.Dekking investigated the limit properties of curves generated by L-systems [32] and concentrated on the problem of determining the fractal(Hausdorff) dimension of the limit set [31]. Prusinkiewicz focused onan interpretation based on a LOGO-style turtle [1] and presented moreexamples of fractals and plant-like structures modeled using L-systems[109, 111]. Further applications of L-systems with turtle interpretationinclude realistic modeling of herbaceous plants [117], description of ko-lam patterns (an art form from Southern India) [112, 115, 133, 134],synthesis of musical scores [110] and automatic generation of space-filling curves [116].

The basic idea of turtle interpretation is given below. A state of theTurtleturtle is defined as a triplet (x, y, α), where the Cartesian coordinates(x, y) represent the turtle’s position, and the angle α, called the heading,is interpreted as the direction in which the turtle is facing. Giventhe step size d and the angle increment δ, the turtle can respond to

1.3. Turtle interpretation of strings 7

Figure 1.5: (a) Turtle interpretation of string symbols F , +, −. (b) Inter-pretation of a string. The angle increment δ is equal to 90◦. Initially theturtle faces up.

commands represented by the following symbols (Figure 1.5a):

F Move forward a step of length d. The state of the turtlechanges to (x′, y′, α), where x′ = x + d cos α and y′ =y + d sin α. A line segment between points (x, y) and(x′, y′) is drawn.

f Move forward a step of length d without drawing a line.

+ Turn left by angle δ. The next state of the turtle is(x, y, α+δ). The positive orientation of angles is counter-clockwise.

− Turn right by angle δ. The next state of the turtle is(x, y, α − δ).

Given a string ν, the initial state of the turtle (x0, y0, α0) and fixed Interpretationparameters d and δ, the turtle interpretation of ν is the figure (set oflines) drawn by the turtle in response to the string ν (Figure 1.5b).Specifically, this method can be applied to interpret strings which aregenerated by L-systems. For example, Figure 1.6 presents four approxi-mations of the quadratic Koch island taken from Mandelbrot’s book [95,page 51]. These figures were obtained by interpreting strings generatedby the following L-system:

ω : F − F − F − Fp : F → F − F + F + FF − F − F + F

The images correspond to the strings obtained in derivations of length0 to 3. The angle increment δ is equal to 90◦. The step length d isdecreased four times between subsequent images, making the distance

8 Chapter 1. Graphical modeling using L-systems

Figure 1.6: Generating a quadratic Koch island

between the endpoints of the successor polygon equal to the length ofthe predecessor segment.

The above example reveals a close relationship between Koch con- Kochconstructionsvs. L-systems

structions and L-systems. The initiator corresponds to the axiom andthe generator corresponds to the production successor. The predeces-sor F represents a single edge. L-systems specified in this way can beperceived as codings for Koch constructions. Figure 1.7 presents furtherexamples of Koch curves generated using L-systems. A slight compli-cation occurs if the curve is not connected; a second production (withthe predecessor f) is then required to keep components the proper dis-tance from each other (Figure 1.8). The ease of modifying L-systemsmakes them suitable for developing new Koch curves. For example, onecan start from a particular L-system and observe the results of insert-ing, deleting or replacing some symbols. A variety of curves obtainedthis way are shown in Figure 1.9.

1.3. Turtle interpretation of strings 9

Figure 1.7: Examples of Koch curves generated using L-systems: (a)Quadratic Koch island [95, page 52], (b) A quadratic modification of thesnowflake curve [95, page 139]

Figure 1.8: Combination of islands and lakes [95, page 121]

10 Chapter 1. Graphical modeling using L-systems

Figure 1.9: A sequence of Koch curves obtained by successive modificationof the production successor

1.4. Synthesis of DOL-systems 11

a n=10, δ=90◦FlFl→Fl+Fr+Fr→-Fl-Fr

b n=6, δ=60◦FrFl→Fr+Fl+FrFr→Fl-Fr-Fl

Figure 1.10: Examples of curves generated by edge-rewriting L-systems: (a)the dragon curve [48], (b) the Sierpinski gasket [132]

1.4 Synthesis of DOL-systems

Random modification of productions gives little insight into the rela-tionship between L-systems and the figures they generate. However,we often wish to construct an L-system which captures a given struc-ture or sequence of structures representing a developmental process.This is called the inference problem in the theory of L-systems. Al-though some algorithms for solving it were reported in the literature[79, 88, 89], they are still too limited to be of practical value in themodeling of higher plants. Consequently, the methods introduced be-low are more intuitive in nature. They exploit two modes of operationfor L-systems with turtle interpretation, called edge rewriting and noderewriting using terminology borrowed from graph grammars [56, 57, 87].In the case of edge rewriting, productions substitute figures for poly-gon edges, while in node rewriting, productions operate on polygonvertices. Both approaches rely on capturing the recursive structure offigures and relating it to a tiling of a plane. Although the concepts areillustrated using abstract curves, they apply to branching structuresfound in plants as well.

1.4.1 Edge rewriting

Edge rewriting can be viewed as an extension of Koch constructions.For example, Figure 1.10a shows the dragon curve [21, 48, 95] and theL-system that generated it. Both the Fl and Fr symbols representedges created by the turtle executing the “move forward” command.The productions substitute Fl or Fr edges by pairs of lines forming

12 Chapter 1. Graphical modeling using L-systems

an=4, δ=60◦FlFl→Fl+Fr++Fr-Fl--FlFl-Fr+Fr→-Fl+FrFr++Fr+Fl--Fl-Fr

bn=2, δ=90◦-FrFl→FlFl-Fr-Fr+Fl+Fl-Fr-FrFl+

Fr+FlFlFr-Fl+Fr+FlFl+Fr-FlFr-Fr-Fl+Fl+FrFr-

Fr→+FlFl-Fr-Fr+Fl+FlFr+Fl-FrFr-Fl-Fr+FlFrFr-Fl-FrFl+Fl+Fr-Fr-Fl+Fl+FrFr

Figure 1.11: Examples of FASS curves generated by edge-rewriting L-systems: (a) hexagonal Gosper curve [51], (b) quadratic Gosper curve [32]or E-curve [96]

left or right turns. Many interesting curves can be obtained assumingtwo types of edges, “left” and “right.” Figures 1.10b and 1.11 presentadditional examples.

The curves included in Figure 1.11 belong to the class of FASSFASS curveconstruction curves (an acronym for space-filling, self-avoiding, simple and self-

similar) [116], which can be thought of as finite, self-avoiding approxi-mations of curves that pass through all points of a square (space-fillingcurves [106]). McKenna [96] presented an algorithm for constructingFASS curves using edge replacement. It exploits the relationship be-tween such a curve and a recursive subdivision of a square into tiles.For example, Figure 1.12 shows the tiling that corresponds to the E-curve of Figure 1.11b. The polygon replacing an edge Fl (Figure 1.12a)approximately fills the square on the left side of Fl (b). Similarly, thepolygon replacing an edge Fr (c) approximately fills the square on theright side of that edge (d). Consequently, in the next derivation step,each of the 25 tiles associated with the curves (b) or (d) will be coveredby their reduced copies (Figure 1.11b). A recursive application of thisargument indicates that the whole curve is approximately space-filling.It is also self-avoiding due to the following two properties:

1.4. Synthesis of DOL-systems 13

Fl→FlFl+Fr+Fr-Fl-Fl+Fr+FrFl-Fr-FlFlFr+Fl-Fr-FlFl-Fr+FlFr+Fr+Fl-Fl-FrFr+

Fr→-FlFl+Fr+Fr-Fl-FlFr-Fl+FrFr+Fl+Fr-FlFrFr+Fl+FrFl-Fl-Fr+Fr+Fl-Fl-FrFr

Figure 1.12: Construction of the E-curve on the square grid. Left and rightedges are distinguished by the direction of ticks.

• the generating polygon is self-avoiding, and

• no matter what the relative orientation of the polygons lying ontwo adjacent tiles, their union is a self-avoiding curve.

The first property is obvious, while the second can be verified by con-sidering all possible relative positions of a pair of adjacent tiles.

Using a computer program to search the space of generating poly-gons, McKenna found that the E-curve is the simplest FASS curveobtained by edge replacement in a square grid. Other curves requiregenerators with more edges (Figure 1.13). The relationship betweenedge rewriting and tiling of the plane extends to branching structures,providing a method for constructing and analyzing L-systems whichoperate according to the edge-rewriting paradigm (see Section 1.10.3).

1.4.2 Node rewriting

The idea of node rewriting is to substitute new polygons for nodes of the Subfigurespredecessor curve. In order to make this possible, turtle interpretationis extended by symbols which represent arbitrary subfigures. As shownin Figure 1.14, each subfigure A from a set of subfigures A is representedby:

• two contact points, called the entry point PA and the exit pointQA, and

• two direction vectors, called the entry vector �pA and the exit vector�qA.

During turtle interpretation of a string ν, a symbol A ∈ A incorporatesthe corresponding subfigure into a picture. To this end, A is translated

14 Chapter 1. Graphical modeling using L-systems

Figure 1.13: Examples of FASS curves generated on the square grid usingedge replacement: (a) a SquaRecurve (grid size 7 × 7), (b) an E-tour (gridsize 9 × 9). Both curves are from [96].

Figure 1.14: Description of a subfigure A

1.4. Synthesis of DOL-systems 15

Ln RnLn+1 Rn+1

Ln+2 Rn+2

Figure 1.15: Recursive construction of the Hilbert curve [63] in terms ofnode replacement

and rotated in order to align its entry point PA and direction �pA withthe current position and orientation of the turtle. Having placed A, theturtle is assigned the resulting exit position QA and direction �qA.

For example, assuming that the contact points and directions of Recursiveformulassubfigures Ln and Rn are as in Figure 1.15, the figures Ln+1 and Rn+1

are captured by the following formulas:

Ln+1 = +RnF − LnFLn − FRn+Rn+1 = −LnF + RnFRn + FLn−

Suppose that curves L0 and R0 are given. One way of evaluating thestring Ln (or Rn) for n > 0 is to generate successive strings recur-sively, in the order of decreasing value of index n. For example, thecomputation of L2 would proceed as follows:

L2 = +R1F − L1FL1 − FR1+= +(−L0F + R0FR0 + FL0−)F − (+R0F − L0FL0 − FR0+)

F (+R0F − L0FL0 − FR0+) − F (−L0F + R0FR0 + FL0−)+

16 Chapter 1. Graphical modeling using L-systems

Thus, the generation of string Ln can be considered as a string-rewritingmechanism, where the symbols on the left side of the recursive formulasare substituted by corresponding strings on the right side. The substi-tution proceeds in a parallel way with F, + and − replacing themselves.Since all indices in any given string have the same value, they can bedropped, provided that a global count of derivation steps is kept. Con-sequently, string Ln can be obtained in a derivation of length n usingthe following L-system:

ω : Lp1 : L → +RF − LFL − FR+p2 : R → −LF + RFR + FL−

In order to complete the curve definition, it is necessary to specifythe subfigures represented by symbols L and R. In the special case ofPure curvespure curves [116], these subfigures are reduced to single points. Thus,one can assume that symbols L and R are erased (replaced by the emptystring) at the end of the derivation. Alternatively, they can be left inthe string and ignored by the turtle during string interpretation. Thissecond approach is consistent with previous definitions of turtle inter-pretation [109, 112]. A general discussion of the relationship betweenrecurrent formulas and L-systems is presented in [61, 62].

Construction of the L-system generating the Hilbert curve can beFASS curveconstruction extended to other FASS curves [116]. Consider an array of m×m square

tiles, each including a smaller square, called a frame. The edges of theframe run at some distance from the tile’s edges. Each frame bounds anopen self-avoiding polygon. The endpoints of this polygon coincide withthe two contact vertices of the frame. Suppose that a single-stroke linerunning through all tiles can be constructed by connecting the contactvertices of neighboring frames using short horizontal or vertical linesegments. A FASS curve can be constructed by the recursive repetitionof this connecting pattern. To this end, the array of m × m connectedtiles is considered a macrotile which contains an open polygon inscribedinto a macroframe. An array of m × m macrotiles is formed, and thepolygons inscribed into the macroframes are connected together. Thisconstruction is carried out recursively, with m × m macrotiles at leveln yielding one macrotile at level n + 1.

Tile arrangements suitable for the generation of FASS curves canbe found algorithmically, by searching the space of all possible arrange-ments on a grid of a given size. Examples of curves synthesized thisway are given in Figures 1.16 and 1.17.

As in the case of edge rewriting, the relationship between noderewriting and tilings of the plane extends to branching structures. Itoffers a method for synthesizing L-systems that generate objects with agiven recursive structure, and links methods for plant generation basedon L-systems with those using iterated function systems [7] (see Chap-ter 8).

1.4. Synthesis of DOL-systems 17

an=3,δ=90◦-LL→LF+RFR+FL-F-LFLFL-FRFR+R→-LFLF+RFRFR+F+RF-LFL-FR

bn=2, δ=90◦-LL→LFLF+RFR+FLFL-FRF-LFL-

FR+F+RF-LFL-FRFRFR+R→-LFLFLF+RFR+FL-F-LF+RFR+

FLF+RFRF-LFL-FRFR

Figure 1.16: Sample FASS curves constructed using tiles with contact pointspositioned along a tile edge: (a) 3 × 3 tiles form a macrotile, (b) 4 × 4 tilesform a macrotile

a n=2, δ=90◦LL→LFRFL-F-RFLFR+F+LFRFLR→RFLFR+F+LFRFL-F-RFLFR

b n=2, δ=45◦LL→L+F+R-F-L+F+R-F-L-F-R+F+L-F-R-F-L+F+R-F-L-F-R-F-

L+F+R+F+L+F+R-F-L+F+R+F+L-R-F+F+L+F+R-F-L+F+R-F-LR→R-F-L+F+R-F-L+F+R+F+L-F-R+F+L+F+R-F-L+F+R+F+L+F+

R-F-L-F-R-F-L+F+R-F-L-F-R+F+L-F-R-F-L+F+R-F-L+F+R

Figure 1.17: Sample FASS curves constructed using tiles with contact pointspositioned diagonally: (a) 3 × 3 tiles form a macrotile (Peano curve [106]),(b) 5 × 5 tiles form a macrotile

18 Chapter 1. Graphical modeling using L-systems

1.4.3 Relationship between edge andnode rewriting

The classes of curves that can be generated using the edge-rewriting andnode-rewriting techniques are not disjoint. For example, reconsider theL-system that generates the dragon curve using edge replacement:

ω : Fl

p1 : Fl → Fl + Fr+p2 : Fr → −Fl − Fr

Assume temporarily that a production predecessor can contain morethan one letter; thus an entire subword can be replaced by the successorof a single production (a formalization of this concept is termed pseudo-L-systems [109]). The dragon-generating L-system can be rewritten as:

ω : Flp1 : Fl → Fl + rF+p2 : rF → −Fl − rF

where the symbols l and r are not interpreted by the turtle. Productionp1 replaces the letter l by the string l + rF− while the leading letter Fis left intact. In a similar way, production p2 replaces the letter r bythe string −Fl−r and leaves the trailing F intact. Thus, the L-systemcan be transformed into node-rewriting form as follows:

ω : Flp1 : l → l + rF+p2 : r → −Fl − r

In practice, the choice between edge rewriting and node rewritingis often a matter of convenience. Neither approach offers an auto-matic, general method for constructing L-systems that capture givenstructures. However, the distinction between edge and node rewritingmakes it easier to understand the intricacies of L-system operation, andin this sense assists in the modeling task. Specifically, the problem offilling a region by a self-avoiding curve is biologically relevant, sincesome plant structures, such as leaves, may tend to fill a plane withoutoverlapping [38, 66, 67, 94].

1.5 Modeling in three dimensions

Turtle interpretation of L-systems can be extended to three dimensionsfollowing the ideas of Abelson and diSessa [1]. The key concept is torepresent the current orientation of the turtle in space by three vectors�H,�L, �U, indicating the turtle’s heading, the direction to the left, and thedirection up. These vectors have unit length, are perpendicular to each

1.5. Modeling in three dimensions 19

Figure 1.18: Controlling the turtle in three dimensions

other, and satisfy the equation �H × �L = �U. Rotations of the turtle arethen expressed by the equation

[�H ′ �L′ �U ′

]=

[�H �L �U

]R,

where R is a 3× 3 rotation matrix [40]. Specifically, rotations by angleα about vectors �U,�L and �H are represented by the matrices:

RU(α) =

cos α sin α 0− sin α cos α 0

0 0 1

RL(α) =

cos α 0 − sin α

0 1 0sin α 0 cos α

RH(α) =

1 0 0

0 cos α − sin α0 sin α cos α

The following symbols control turtle orientation in space (Figure 1.18):

+ Turn left by angle δ, using rotation matrix RU(δ).

− Turn right by angle δ, using rotation matrix RU(−δ).

& Pitch down by angle δ, using rotation matrix RL(δ).

∧ Pitch up by angle δ, using rotation matrix RL(−δ).

\ Roll left by angle δ, using rotation matrix RH(δ).

/ Roll right by angle δ, using rotation matrix RH(−δ).

| Turn around, using rotation matrix RU(180◦).

20 Chapter 1. Graphical modeling using L-systems

n=2, δ=90◦AA → B-F+CFC+F-D&F∧D-F+&&CFC+F+B//B → A&F∧CFB∧F∧D∧∧-F-D∧|F∧B|FC∧F∧A//C → |D∧|F∧B-F+C∧F∧A&&FA&F∧C+F+B∧F∧D//D → |CFB-F+B|FA&F∧A&&FB-F+B|FC//

Figure 1.19: A three-dimensional extension of the Hilbert curve [139]. Col-ors represent three-dimensional “frames” associated with symbols A (red), B(blue), C (green) and D (yellow).

1.6. Branching structures 21

As an example of a three-dimensional object created using an L-system, consider the extension of the Hilbert curve shown in Figure 1.19.The L-system was constructed with the node-replacement techniquediscussed in the previous section, using cubes and “macrocubes” in-stead of tiles and macrotiles.

1.6 Branching structures

According to the rules presented so far, the turtle interprets a characterstring as a sequence of line segments. Depending on the segment lengthsand the angles between them, the resulting line is self-intersecting ornot, can be more or less convoluted, and may have some segmentsdrawn many times and others made invisible, but it always remainsjust a single line. However, the plant kingdom is dominated by branch-ing structures; thus a mathematical description of tree-like shapes andmethods for generating them are needed for modeling purposes. Anaxial tree [89, 117] complements the graph-theoretic notion of a rootedtree [108] with the botanically motivated notion of branch axis.

1.6.1 Axial trees

A rooted tree has edges that are labeled and directed. The edge se-quences form paths from a distinguished node, called the root or base,to the terminal nodes. In the biological context, these edges are re-ferred to as branch segments . A segment followed by at least one moresegment in some path is called an internode. A terminal segment (withno succeeding edges) is called an apex.

An axial tree is a special type of rooted tree (Figure 1.20). At eachof its nodes, at most one outgoing straight segment is distinguished.All remaining edges are called lateral or side segments. A sequence ofsegments is called an axis if:

• the first segment in the sequence originates at the root of the treeor as a lateral segment at some node,

• each subsequent segment is a straight segment, and

• the last segment is not followed by any straight segment in thetree.

Together with all its descendants, an axis constitutes a branch. Abranch is itself an axial (sub)tree.

Axes and branches are ordered. The axis originating at the root ofthe entire plant has order zero. An axis originating as a lateral segmentof an n-order parent axis has order n+1. The order of a branch is equalto the order of its lowest-order or main axis.

22 Chapter 1. Graphical modeling using L-systems

Figure 1.20: An axial tree

Figure 1.21: Sample tree generated using a method based on Horton–Strahler analysis of branching patterns

1.6. Branching structures 23

Figure 1.22: A tree production p and its application to the edge S in a treeT1

Axial trees are purely topological objects. The geometric connotationof such terms as straight segment, lateral segment and axis should beviewed at this point as an intuitive link between the graph-theoreticformalism and real plant structures.

The proposed scheme for ordering branches in axial trees was in-troduced originally by Gravelius [53]. MacDonald [94, pages 110–121]surveys this and other methods applicable to biological and geograph-ical data such as stream networks. Of special interest are methodsproposed by Horton [70, 71] and Strahler, which served as a basis forsynthesizing botanical trees [37, 152] (Figure 1.21).

1.6.2 Tree OL-systems

In order to model development of branching structures, a rewritingmechanism can be used that operates directly on axial trees. A rewrit-ing rule, or tree production, replaces a predecessor edge by a successoraxial tree in such a way that the starting node of the predecessor isidentified with the successor’s base and the ending node is identifiedwith the successor’s top (Figure 1.22).

A tree OL-system G is specified by three components: a set of edgelabels V , an initial tree ω with labels from V , and a set of tree produc-tions P . Given the L-system G, an axial tree T2 is directly derived froma tree T1, noted T1 ⇒ T2, if T2 is obtained from T1 by simultaneouslyreplacing each edge in T1 by its successor according to the productionset P . A tree T is generated by G in a derivation of length n if thereexists a sequence of trees T0, T1, . . . , Tn such that T0 = ω, Tn = T andT0 ⇒ T1 ⇒ . . . ⇒ Tn.

24 Chapter 1. Graphical modeling using L-systems

Figure 1.23: Bracketed string representation of an axial tree

1.6.3 Bracketed OL-systems

The definition of tree L-systems does not specify the data structure forrepresenting axial trees. One possibility is to use a list representationwith a tree topology. Alternatively, axial trees can be represented usingstrings with brackets [82]. A similar distinction can be observed in Kochconstructions, which can be implemented either by rewriting edges andpolygons or their string representations. An extension of turtle in-terpretation to strings with brackets and the operation of bracketedL-systems [109, 111] are described below.

Two new symbols are introduced to delimit a branch. They areinterpreted by the turtle as follows:

[ Push the current state of the turtle onto a pushdownStackoperations stack. The information saved on the stack contains the

turtle’s position and orientation, and possibly other at-tributes such as the color and width of lines being drawn.

] Pop a state from the stack and make it the current stateof the turtle. No line is drawn, although in general theposition of the turtle changes.

An example of an axial tree and its string representation are shownin Figure 1.23.

Derivations in bracketed OL-systems proceed as in OL-systems with-2D structuresout brackets. The brackets replace themselves. Examples of two-dimensional branching structures generated by bracketed OL-systemsare shown in Figure 1.24.

Figure 1.25 is an example of a three-dimensional bush-like structureBush-likestructure generated by a bracketed L-system. Production p1 creates three new

branches from an apex of the old branch. A branch consists of anedge F forming the initial internode, a leaf L and an apex A (whichwill subsequently create three new branches). Productions p2 and p3

1.6. Branching structures 25

an=5,δ=25.7◦FF→F[+F]F[-F]F

bn=5,δ=20◦FF→F[+F]F[-F][F]

cn=4,δ=22.5◦FF→FF-[-F+F+F]+

[+F-F-F]

dn=7,δ=20◦XX→F[+X]F[-X]+XF→FF

en=7,δ=25.7◦XX→F[+X][-X]FXF→FF

fn=5,δ=22.5◦XX→F-[[X]+X]+F[+FX]-XF→FF

Figure 1.24: Examples of plant-like structures generated by bracketed OL-systems. L-systems (a), (b) and (c) are edge-rewriting, while (d), (e) and(f) are node-rewriting.

26 Chapter 1. Graphical modeling using L-systems

n=7, δ=22.5◦

ω : Ap1 : A → [&FL!A]/////’[&FL!A]///////’[&FL!A]p2 : F → S ///// Fp3 : S → F Lp4 : L → [’’’∧∧{-f+f+f-|-f+f+f}]

Figure 1.25: A three-dimensional bush-like structure

specify internode growth. In subsequent derivation steps the internodegets longer and acquires new leaves. This violates a biological ruleof subapical growth (discussed in detail in Chapter 3), but producesan acceptable visual effect in a still picture. Production p4 specifiesthe leaf as a filled polygon with six edges. Its boundary is formedfrom the edges f enclosed between the braces { and } (see Chapter 5for further discussion). The symbols ! and ′ are used to decrementthe diameter of segments and increment the current index to the colortable, respectively.

Another example of a three-dimensional plant is shown in Fig-Plantwith flowers ure 1.26. The L-system can be described and analyzed in a way similar

to the previous one.

1.6. Branching structures 27

n=5, δ=18◦

ω : plantp1 : plant → internode + [ plant + flower] − − //

[ − − leaf ] internode [ + + leaf ] −[ plant flower ] + + plant flower

p2 : internode → F seg [// & & leaf ] [// ∧ ∧ leaf ] F segp3 : seg → seg F segp4 : leaf → [’ { +f−ff−f+ | +f−ff−f } ]p5 : flower → [ & & & pedicel ‘ / wedge //// wedge ////

wedge //// wedge //// wedge ]p6 : pedicel → FFp7 : wedge → [‘ ∧ F ] [ { & & & & −f+f | −f+f } ]

Figure 1.26: A plant generated by an L-system

28 Chapter 1. Graphical modeling using L-systems

1.7 Stochastic L-systems

All plants generated by the same deterministic L-system are identical.An attempt to combine them in the same picture would produce astriking, artificial regularity. In order to prevent this effect, it is nec-essary to introduce specimen-to-specimen variations that will preservethe general aspects of a plant but will modify its details.

Variation can be achieved by randomizing the turtle interpretation,the L-system, or both. Randomization of the interpretation alone hasa limited effect. While the geometric aspects of a plant — such asthe stem lengths and branching angles — are modified, the underly-ing topology remains unchanged. In contrast, stochastic applicationof productions may affect both the topology and the geometry of theplant. The following definition is similar to that of Yokomori [162] andEichhorst and Savitch [35].

A stochastic OL-system is an ordered quadruplet Gπ = 〈V, ω, P, π〉.L-systemThe alphabet V , the axiom ω and the set of productions P are definedas in an OL-system (page 4). Function π : P → (0, 1], called theprobability distribution, maps the set of productions into the set ofproduction probabilities. It is assumed that for any letter a ∈ V , thesum of probabilities of all productions with the predecessor a is equalto 1.

We will call the derivation µ ⇒ ν a stochastic derivation in Gπ if forDerivationeach occurrence of the letter a in the word µ the probability of applyingproduction p with the predecessor a is equal to π(p). Thus, differentproductions with the same predecessor can be applied to various occur-rences of the same letter in one derivation step.

A simple example of a stochastic L-system is given below.Example

ω : F

p1 : F.33→ F [+F ]F [−F ]F

p2 : F.33→ F [+F ]F

p3 : F.34→ F [−F ]F

The production probabilities are listed above the derivation symbol→. Each production can be selected with approximately the sameprobability of 1/3. Examples of branching structures generated by thisL-system with derivations of length 5 are shown in Figure 1.27. Notethat these structures look like different specimens of the same (albeitfictitious) plant species.

A more complex example is shown in Figure 1.28. The field consistsFlower fieldof four rows and four columns of plants. All plants are generated by astochastic modification of the L-system used to generate Figure 1.26.

1.7. Stochastic L-systems 29

Figure 1.27: Stochastic branching structures

Figure 1.28: Flower field

30 Chapter 1. Graphical modeling using L-systems

The essence of this modification is to replace the original productionp3 by the following three productions:

p′3 : seg.33→ seg [ // & & leaf ] [// ∧∧ leaf ] F seg

p′′3 : seg.33→ seg F seg

p′′′3 : seg.34→ seg

Thus, in any step of the derivation, the stem segment seg may eithergrow and produce new leaves (production p′3), grow without producingnew leaves (production p′′3), or not grow at all (production p′′′3 ). All threeevents occur with approximately the same probability. The resultingfield appears to consist of various specimens of the same plant species.If the same L-system was used again (with different seed values for therandom number generator), a variation of this image would be obtained.

1.8 Context-sensitive L-systems

Productions in OL-systems are context-free; i.e. applicable regardlessContext instringL-systems

of the context in which the predecessor appears. However, productionapplication may also depend on the predecessor’s context. This effect isuseful in simulating interactions between plant parts, due for example tothe flow of nutrients or hormones. Various context-sensitive extensionsof L-systems have been proposed and studied thoroughly in the past[62, 90, 128]. 2L-systems use productions of the form al < a > ar → χ,where the letter a (called the strict predecessor) can produce word χ ifand only if a is preceded by letter al and followed by ar. Thus, lettersal and ar form the left and the right context of a in this production.Productions in 1L-systems have one-sided context only; consequently,they are either of the form al < a → χ or a > ar → χ. OL-systems,1L-systems and 2L-systems belong to a wider class of IL-systems, alsocalled (k,l)-systems. In a (k,l)-system, the left context is a word oflength k and the right context is a word of length l.

In order to keep specifications of L-systems short, the usual notionof IL-systems has been modified here by allowing productions withdifferent context lengths to coexist within a single system. Further-more, context-sensitive productions are assumed to have precedenceover context-free productions with the same strict predecessor. Conse-quently, if a context-free and a context-sensitive production both applyto a given letter, the context-sensitive one should be selected. If no pro-duction applies, this letter is replaced by itself as previously assumedfor OL-systems.

1.8. Context-sensitive L-systems 31

Figure 1.29: The predecessor of a context-sensitive tree production (a)matches edge S in a tree T (b)

The following sample 1L-system makes use of context to simulate signal Signalpropagationpropagation throughout a string of symbols:

ω : baaaaaaaap1 : b < a → bp2 : b → a

The first few words generated by this L-system are given below:

baaaaaaaaabaaaaaaaaabaaaaaaaaabaaaaaaaaabaaaa· · ·

The letter b moves from the left side to the right side of the string.A context-sensitive extension of tree L-systems requires neighbor Context in tree

L-systemsedges of the replaced edge to be tested for context matching. A prede-cessor of a context-sensitive production p consists of three components:a path l forming the left context, an edge S called the strict predecessor,and an axial tree r constituting the right context (Figure 1.29). Theasymmetry between the left context and the right context reflects thefact that there is only one path from the root of a tree to a given edge,while there can be many paths from this edge to various terminal nodes.Production p matches a given occurrence of the edge S in a tree T if lis a path in T terminating at the starting node of S, and r is a subtree

32 Chapter 1. Graphical modeling using L-systems

of T originating at the ending node of S. The production can then beapplied by replacing S with the axial tree specified as the productionsuccessor.

The introduction of context to bracketed L-systems is more difficultContext inbracketedL-systems

than in L-systems without brackets, because the bracketed string repre-sentation of axial trees does not preserve segment neighborhood. Conse-quently, the context matching procedure may need to skip over symbolsrepresenting branches or branch portions. For example, Figure 1.29 in-dicates that a production with the predecessor BC < S > G[H]M canbe applied to symbol S in the string

ABC[DE][SG[HI[JK]L]MNO],

which involves skipping over symbols [DE] in the search for left context,and I[JK]L in the search for right context.

Within the formalism of bracketed L-systems, the left context canbe used to simulate control signals that propagate acropetally, i.e., fromthe root or basal leaves towards the apices of the modeled plant, whilethe right context represents signals that propagate basipetally, i.e., fromthe apices towards the root. For example, the following 1L-systemsimulates propagation of an acropetal signal in a branching structurethat does not grow:

#ignore : +−

ω : Fb[+Fa]Fa[−Fa]Fa[+Fa]Fa

p1 : Fb < Fa → Fb

Symbol Fb represents a segment already reached by the signal, whileFa represents a segment that has not yet been reached. The #ignorestatement indicates that the geometric symbols + and − should beconsidered as non-existent while context matching. Images representingconsecutive stages of signal propagation (corresponding to consecutivewords generated by the L-system under consideration) are shown inFigure 1.30a.

The propagation of a basipetal signal can be simulated in a similarway (Figure 1.30b):

#ignore : +−

ω : Fa[+Fa]Fa[−Fa]Fa[+Fa]Fb

p1 : Fa > Fb → Fb

1.8. Context-sensitive L-systems 33

Figure 1.30: Signal propagation in a branching structure: (a) acropetal, (b)basipetal

The operation of context-sensitive L-systems is examined further using L-systems ofHogeweg,Hesper andSmith

examples obtained by Hogeweg and Hesper [64]. In 1974, they pub-lished the results of an exhaustive study of 3,584 patterns generatedby a class of bracketed 2L-systems defined over the alphabet {0,1}.Some of these patterns had plant-like shapes. Subsequently, Smithsignificantly improved the quality of the generated images using state-of-the-art computer imagery techniques [136, 137]. Sample structuresgenerated by L-systems similar to those proposed by Hogeweg and Hes-per are shown in Figure 1.31. The differences are related to the geo-metric interpretation of the resulting strings. According to the originalinterpretation, consecutive branches are issued alternately to the leftand right, whereas turtle interpretation requires explicit specificationof branching angles within the L-system.

34 Chapter 1. Graphical modeling using L-systems

Figure 1.31: Examples of branching structures generated using L-systemsbased on the results of Hogeweg and Hesper [64]

1.8. Context-sensitive L-systems 35

a n=30,δ=22.5◦#ignore: +-FF1F1F10 < 0 > 0 → 00 < 0 > 1 → 1[+F1F1]0 < 1 > 0 → 10 < 1 > 1 → 11 < 0 > 0 → 01 < 0 > 1 → 1F11 < 1 > 0 → 01 < 1 > 1 → 0* < + > * → -* < - > * → +

b n=30,δ=22.5◦#ignore: +-FF1F1F10 < 0 > 0 → 10 < 0 > 1 → 1[-F1F1]0 < 1 > 0 → 10 < 1 > 1 → 11 < 0 > 0 → 01 < 0 > 1 → 1F11 < 1 > 0 → 11 < 1 > 1 → 0* < + > * → -* < - > * → +

c n=26, δ=25.75◦#ignore: +-FF1F1F10 < 0 > 0 → 00 < 0 > 1 → 10 < 1 > 0 → 00 < 1 > 1 → 1[+F1F1]1 < 0 > 0 → 01 < 0 > 1 → 1F11 < 1 > 0 → 01 < 1 > 1 → 0* < - > * → +* < + > * → -

d n=24, δ=25.75◦#ignore: +-FF0F1F10 < 0 > 0 → 10 < 0 > 1 → 00 < 1 > 0 → 00 < 1 > 1 → 1F11 < 0 > 0 → 11 < 0 > 1 → 1[+F1F1]1 < 1 > 0 → 11 < 1 > 1 → 0* < + > * → -* < - > * → +

e n=26, δ=22.5◦#ignore: +-FF1F1F10 < 0 > 0 → 00 < 0 > 1 → 1[-F1F1]0 < 1 > 0 → 10 < 1 > 1 → 11 < 0 > 0 → 01 < 0 > 1 → 1F11 < 1 > 0 → 11 < 1 > 1 → 0* < + > * → -* < - > * → +

Figure 1.31 (continued): L-systems of Hogeweg and Hesper

36 Chapter 1. Graphical modeling using L-systems

1.9 Growth functions

During the synthesis of a plant model it is often convenient to dis-Exponentialgrowth tinguish productions that specify the branching pattern from those

that describe elongation of plant segments. This separation can beobserved in some of the L-systems considered so far. For example, inL-systems (d), (e) and (f) from Figure 1.24 the first productions cap-ture the branching patterns, while the remaining productions, equal inall cases to F → FF , describe elongation of segments represented bysequences of symbols F . The number of letters F in a string χn orig-inating from a single letter F is doubled in each derivation step, thusthe elongation is exponential, with length(χn) = 2n.

A function that describes the number of symbols in a word in termsBasicproperties of its derivation length is called a growth function. The theory of L-

systems contains an extensive body of results on growth functions [62,127]. The central observation is that the growth functions of DOL-systems are independent of the letter ordering in the productions andderived words. Consequently, the relation between the number of letteroccurrences in a pair of words µ and ν, such that µ ⇒ ν, can beconveniently expressed using matrix notation.

Let G = 〈V, ω, P 〉 be a DOL-system and assume that letters ofthe alphabet V have been ordered, V = {a1, a2, . . . , am}. Construct asquare matrix Qm×m, where entry qij is equal to the number of occur-rences of letter aj in the successor of the production with predecessorai. Let ak

i denote the number of occurrences of letter ai in the wordx generated by G in a derivation of length k. The definition of directderivation in a DOL-system implies that

[ak

1 ak2 · · · ak

m

]

q11 q12 · · · q1m

q21 q22 · · · q2m...

qm1 qm2 · · · qmm

=

[ak+1

1 ak+12 · · · ak+1

m

].

This matrix notation is useful in the analysis of growth functions. Forexample, consider the following L-system:

ω : ap1 : a → abp2 : b → a

(1.2)

The relationship between the number of occurrences of letters a and bin two consecutively derived words is

[ak bk

] [1 11 0

]=

[ak+1 bk+1

]

1.9. Growth functions 37

orak+1 = ak + bk = ak + ak−1

for k = 1, 2, 3, . . . . From the axiom it follows that a0 = 1 anda1 = b0 = 0. Thus, the number of letters a in the strings gener-ated by the L-system specified in equation (1.2) grows according to theFibonacci series: 1, 1, 2, 3, 5, 8, . . . . This growth function was imple-mented by productions p2 and p3 in the L-system generating the bushin Figure 1.25 (page 26) to describe the elongation of its internodes.

Polynomial growth functions of arbitrary degree can be obtained Polynomialgrowthusing L-systems of the following form:

ω : a0

p1 : a0 → a0a1

p2 : a1 → a1a2

p3 : a2 → a2a3

p4 : a3 → a3a4...

The matrix Q is given below:

Q =

1 1 0 0 · · ·0 1 1 0 · · ·0 0 1 1 · · ·0 0 0 1 · · ·

...

Thus, for any i, k ≥ 1, the number aki of occurrences of symbol ai in

the string generated in a derivation of length k satisfies the equality

aki + ak

i+1 = ak+1i+1 .

Taking into consideration the axiom, the distribution of letters ai asa function of the derivation length is captured by the following table(only non-zero terms are shown):

k ak0 ak

1 ak2 ak

3 ak4 ak

5 ak6 ak

7

0 11 1 12 1 2 13 1 3 3 14 1 4 6 4 15 1 5 10 10 5 16 1 6 15 20 15 6 17 1 7 21 35 35 21 7 1

...

38 Chapter 1. Graphical modeling using L-systems

This table represents the Pascal triangle, thus for any k ≥ i ≥ 1 itsterms satisfy the following equality:

aki =

(ki

)=

k(k − 1) · · · (k − i + 1)

1 · 2 · · · i

Consequently, the number of occurrences of letter ai as a function ofthe derivation length k is expressed by a polynomial of degree i. Byidentifying letter ai with the turtle symbol F , it is possible to model in-ternode elongation expressed by polynomials of arbitrary degree i ≥ 0.This observation was generalized by Szilard [140], who developed an al-gorithm for constructing a DOL-system with growth functions specifiedby any positive, nondecreasing polynomials with integer coefficients [62,page 276].

The examples of growth functions considered so far include expo-Characterizationnential and polynomial functions. Rozenberg and Salomaa [127, pages30–38] show that, in general, the growth function fG(n) of any DOL-system G = 〈V, ω, P 〉 is a combination of polynomial and exponentialfunctions:

fG(n) =s∑

i=1

Pi(n)ρni for n ≥ n0, (1.3)

where Pi(n) denotes a polynomial with integer coefficients, ρi is a non-negative integer, and n0 is the total number of letters in the alphabetV . Unfortunately, many growth processes observed in nature cannot bedescribed by equation (1.3). Two approaches are then possible withinthe framework of the theory of L-systems.

The first is to extend the size n0 of the alphabet V , so that theSigmoidalgrowth growth process of interest will be captured by the initial derivation

steps, ω = µ0 ⇒ µ1 ⇒ · · · ⇒ µn0 , before equation (1.3) starts to apply.For example, the L-system

ω : a0

pi : ai → ai+1b0 for i < kpk+j : bj → bj+1F for j < l

(1.4)

over the alphabet V = {a0, a1, ..., ak}∪{b0, b1, ..., bl}∪{F} can be usedto approximate a sigmoidal elongation of a segment represented by asequence of symbols F (Figure 1.32). The term sigmoidal refers to afunction with a plot in the shape of the letter S. Such functions arecommonly found in biological processes [143], with the initial part ofthe curve representing the growth of a young organism, and the latterpart corresponding to the organism close to its final size.

The second approach to the synthesis of growth functions out-side the class captured by equation (1.3) is to use context-sensitiveSquare-root

growth L-systems. For example, the following 2L-system has a growth func-tion given by fG(n) = �√n + 4, where �x is the floor function.

1.9. Growth functions 39

Figure 1.32: A sigmoidal growth function implemented using the L-systemin equation (1.4), for k = l = 20

ω : XFuFaXp1 : Fu < Fa > Fa → Fu

p2 : Fu < Fa > X → FdFa

p3 : Fa < Fa > Fd → Fd

p4 : X < Fa > Fd → Fu

p5 : Fu → Fa

p6 : Fd → Fa

(1.5)

The operation of this L-system is illustrated in Figure 1.33. Produc-tions p1 and p3, together with p5 and p6, propagate symbols Fu and Fd

up and down the string of symbols µ. Productions p2 and p4 change thepropagation direction, after symbol X marking a string end has beenreached by Fu or Fd, respectively. In addition, p2 extends the stringwith a symbol Fa. Thus, the number of derivation steps increases bytwo between consecutive applications of production p2. As a result,string extension occurs at derivation steps n expressed by the squareof the string length, which yields the growth function �√n + 4.

In practice it is often difficult, if not impossible, to find L-systems Limitationswith the required growth functions. Vitanyi [153] illustrates this byreferring to sigmoidal curves:

If we want to obtain sigmoidal growth curves with theoriginal L-systems then not even the introduction of cellinteraction can help us out. In the first place, we end upconstructing quite unlikely flows of messages through theorganism, which are more suitable to electronic computers,and in fact give the organism universal computing power.Secondly, and this is more fundamental, we can not obtain

40 Chapter 1. Graphical modeling using L-systems

Figure 1.33: Square-root growth implemented using the L-system specifiedin equation (1.5)

growth which, always increasing the size of the organism,tends towards stability in the limit. The slowest increasinggrowth we can obtain by allowing cell interaction is loga-rithmic and thus can not at all account for the asymptoticbehavior of sigmoidal growth functions.

In the next section we present an extension of L-systems that makesit possible to avoid this problem by allowing for explicit inclusion ofgrowth functions into L-system specifications.

1.10 Parametric L-systems

Although L-systems with turtle interpretation make it possible to gen-Motivationerate a variety of interesting objects, from abstract fractals to plant-likebranching structures, their modeling power is quite limited. A majorproblem can be traced to the reduction of all lines to integer multi-ples of the unit segment. As a result, even such a simple figure as anisosceles right-angled triangle cannot be traced exactly, since the ratioof its hypotenuse length to the length of a side is expressed by the irra-tional number

√2. Rational approximation of line length provides only

a limited solution, because the unit step must be the smallest common

��

���

1

1

√2

denominator of all line lengths in the modeled structure. Consequently,the representation of a simple plant module, such as an internode, mayrequire a large number of symbols. The same argument applies to an-gles. Problems become even more pronounced while simulating changesto the modeled structure over time, since some growth functions can-not be expressed conveniently using L-systems. Generally, it is difficult

1.10. Parametric L-systems 41

to capture continuous phenomena, since the obvious technique of dis-cretizing continuous values may require a large number of quantizationlevels, yielding L-systems with hundreds of symbols and productions.Consequently, model specification becomes difficult, and the mathe-matical beauty of L-systems is lost.

In order to solve similar problems, Lindenmayer proposed that nu-merical parameters be associated with L-system symbols [83]. He illus-trated this idea by referring to the continuous development of branchingstructures and diffusion of chemical compounds in a nonbranching fil-ament of Anabaena catenula. Both problems were revisited in laterpapers [25, 77]. A definition of parametric L-systems was formulatedby Prusinkiewicz and Hanan [113] and is presented below.

1.10.1 Parametric OL-systems

Parametric L-systems operate on parametric words, which are strings Parametricwordsof modules consisting of letters with associated parameters. The let-

ters belong to an alphabet V , and the parameters belong to the setof real numbers �. A module with letter A ∈ V and parametersa1, a2, ..., an ∈ � is denoted by A(a1, a2, ..., an). Every module belongsto the set M = V ×�∗, where �∗ is the set of all finite sequences of pa-rameters. The set of all strings of modules and the set of all nonemptystrings are denoted by M∗ = (V ×�∗)∗ and M+ = (V ×�∗)+, respec-tively.

The real-valued actual parameters appearing in the words corre- Expressionsspond with formal parameters used in the specification of L-systemproductions. If Σ is a set of formal parameters, then C(Σ) denotes alogical expression with parameters from Σ, and E(Σ) is an arithmeticexpression with parameters from the same set. Both types of expres-sions consist of formal parameters and numeric constants, combinedusing the arithmetic operators +, −, ∗, /; the exponentiation operator∧, the relational operators <, >, =; the logical operators !, &, | (not,and, or); and parentheses (). Standard rules for constructing syntac-tically correct expressions and for operator precedence are observed.Relational and logical expressions evaluate to zero for false and one fortrue. A logical statement specified as the empty string is assumed tohave value one. The sets of all correctly constructed logical and arith-metic expressions with parameters from Σ are noted C(Σ) and E(Σ).

A parametric OL-system is defined as an ordered quadruplet G = ParametricOL-system〈V, Σ, ω, P 〉, where

• V is the alphabet of the system,

• Σ is the set of formal parameters,

• ω ∈ (V ×�∗)+ is a nonempty parametric word called the axiom,

• P ⊂ (V ×Σ∗)×C(Σ)× (V ×E(Σ))∗ is a finite set of productions.

42 Chapter 1. Graphical modeling using L-systems

The symbols : and → are used to separate the three components of aproduction: the predecessor, the condition and the successor. For exam-ple, a production with predecessor A(t), condition t > 5 and successorB(t + 1)CD(t ∧ 0.5, t − 2) is written as

A(t) : t > 5 → B(t + 1)CD(t ∧ 0.5, t − 2). (1.6)

A production matches a module in a parametric word if the followingDerivationconditions are met:

• the letter in the module and the letter in the production prede-cessor are the same,

• the number of actual parameters in the module is equal to thenumber of formal parameters in the production predecessor, and

• the condition evaluates to true if the actual parameter values aresubstituted for the formal parameters in the production.

A matching production can be applied to the module, creating a stringof modules specified by the production successor. The actual parame-ter values are substituted for the formal parameters according to theirposition. For example, production (1.6) above matches a module A(9),since the letter A in the module is the same as in the production pre-decessor, there is one actual parameter in the module A(9) and oneformal parameter in the predecessor A(t), and the logical expressiont > 5 is true for t = 9. The result of the application of this productionis a parametric word B(10)CD(3, 7).

If a module a produces a parametric word χ as the result of aproduction application in an L-system G, we write a �→ χ. Given aparametric word µ = a1a2...am, we say that the word ν = χ1χ2...χm

is directly derived from (or generated by) µ and write µ =⇒ ν if andonly if ai �→ χi for all i = 1, 2, ...,m. A parametric word ν is generatedby G in a derivation of length n if there exists a sequence of wordsµ0, µ1, ..., µn such that µ0 = ω, µn = ν and µ0 =⇒ µ1 =⇒ ... =⇒ µn.

An example of a parametric L-system is given below.Example

ω : B(2)A(4, 4)p1 : A(x, y) : y <= 3 → A(x ∗ 2, x + y)p2 : A(x, y) : y > 3 → B(x)A(x/y, 0)p3 : B(x) : x < 1 → Cp4 : B(x) : x >= 1 → B(x − 1)

(1.7)

As in the case of non-parametric L-systems, it is assumed that a modulereplaces itself if no matching production is found in the set P . Thewords obtained in the first few derivation steps are shown in Figure 1.34.

1.10. Parametric L-systems 43

Figure 1.34: The initial sequence of strings generated by the parametricL-system specified in equation (1.7)

1.10.2 Parametric 2L-systems

Productions in parametric OL-systems are context-free, i.e., applicableregardless of the context in which the predecessor appears. A context-sensitive extension is necessary to model information exchange betweenneighboring modules. In the parametric case, each component of theproduction predecessor (the left context, the strict predecessor and theright context) is a parametric word with letters from the alphabet Vand formal parameters from the set Σ. Any formal parameters mayappear in the condition and the production successor.

A sample context-sensitive production is given below: Example

A(x) < B(y) > C(z) : x + y + z > 10 → E((x + y)/2)F ((y + z)/2)

It can be applied to the module B(5) that appears in a parametric word

· · ·A(4)B(5)C(6) · · · (1.8)

since the sequence of letters A,B,C in the production and in parametricword (1.8) are the same, the numbers of formal parameters and actualparameters coincide, and the condition 4 + 5 + 6 > 10 is true. As aresult of the production application, the module B(5) will be replacedby a pair of modules E(4.5)F (5.5). Naturally, the modules A(4) andC(6) will be replaced by other productions in the same derivation step.

Parametric 2L-systems provide a convenient tool for expressing de- Anabaena withheterocystsvelopmental models that involve diffusion of substances throughout an

organism. A good example is provided by an extended model of thepattern of cells observed in Anabaena catenula and other blue-greenbacteria [99]. This model was proposed by de Koster and Linden-mayer [25].

44 Chapter 1. Graphical modeling using L-systems

#define CH 900 /* high concentration */#define CT 0.4 /* concentration threshold */#define ST 3.9 /* segment size threshold */#include H /* heterocyst shape specification */#ignore f ∼ H

ω : -(90)F(0,0,CH)F(4,1,CH)F(0,0,CH)

p1 : F(s,t,c) : t=1 & s>=6 →F(s/3*2,2,c)f(1)F(s/3,1,c)

p2 : F(s,t,c) : t=2 & s>=6 →F(s/3,2,c)f(1)F(s/3*2,1,c)

p3 : F(h,i,k) < F(s,t,c) > F(o,p,r) : s>ST|c>CT →F(s+.1,t,c+0.25*(k+r-3*c))

p4 : F(h,i,k) < F(s,t,c) > F(o,p,r) : !(s>ST|c>CT) →F(0,0,CH) ∼ H(1)

p5 : H(s) : s<3 → H(s*1.1)

L-system 1.1: Anabaena catenula

Generally, the bacteria under consideration form a nonbranching fila-ment consisting of two classes of cells: vegetative cells and heterocysts.Usually, the vegetative cells divide and produce two daughter vegeta-tive cells. This mechanism is captured by the L-system specified inequation (1.1) and Figure 1.4 (page 5). However, in some cases thevegetative cells differentiate into heterocysts. Their distribution formsa well-defined pattern, characterized by a relatively constant numberof vegetative cells separating consecutive heterocysts. How does theorganism maintain constant spacing of heterocysts while growing? Themodel explains this phenomenon using a biologically well-motivatedhypothesis that heterocyst distribution is regulated by nitrogen com-pounds produced by the heterocysts, transported from cell to cell acrossthe filament, and decayed in the vegetative cells. If the compound’sconcentration in a young vegetative cell falls below a specific level, thiscell differentiates into a heterocyst (L-system 1.1).

The #define statements assign values to numerical constants usedin the L-system. The #include statement specifies the shape of a het-erocyst (a disk) by referring to a library of predefined shapes (see Sec-tion 5.1). Cells are represented by modules F (s, t, c), where s standsfor cell length, t is cell type (0 - heterocyst, 1 and 2 - vegetative types1),

1The model of Anabaena introduced in Section 1.2 distinguished between fourtypes of cells: ar, br, al and bl. Cells b do not divide and can be considered as youngforms of the corresponding cells a. Thus, the vegetative type 1 considered hereembraces cells ar and br, while type 2 embraces al and bl. The formal relationshipbetween the four-cell and two-cell models is further discussed in Chapter 6.

1.10. Parametric L-systems 45

Figure 1.35: Development of Anabaena catenula with heterocysts, simulatedusing parametric L-system 1.1

and c represents the concentration of nitrogen compounds. Productionsp1 and p2 describe division of the vegetative cells. They each create twodaughter cells of unequal length. The difference between cells of type1 and 2 lies in the ordering of the long and short daughter cells. Pro-duction p3 captures the processes of transportation and decay of thenitrogen compounds. Their concentration c is related to the concen-tration in the neighboring cells k and r, and changes in each derivationstep according to the formula

c′ = c + 0.25(k + r − 3 ∗ c).

Production p4 describes differentiation of a vegetative cell into a hete-rocyst. The condition specifies that this process occurs when the celllength does not exceed the threshold value ST = 3.9 (which meansthat the cell is young enough), and the concentration of the nitrogencompounds falls below the threshold value CT = 0.4. Production p5

describes the subsequent growth of the heterocyst.Snapshots of the developmental sequence of Anabaena are given

in Figure 1.35. The vegetative cells are shown as rectangles, coloredaccording to the concentration of the nitrogen compounds (white meanslow concentration). The heterocysts are represented as red disks. Thevalues of parameters CH, CT and ST were selected to provide the

46 Chapter 1. Graphical modeling using L-systems

correct distribution of the heterocysts, and correspond closely to thevalues reported in [25]. The mathematical model made it possible toestimate these parameters, although they are not directly observable.

1.10.3 Turtle interpretation of parametric words

If one or more parameters are associated with a symbol interpreted bythe turtle, the value of the first parameter controls the turtle’s state.If the symbol is not followed by any parameter, default values specifiedoutside the L-system are used as in the non-parametric case. The basicset of symbols affected by the introduction of parameters is listed below.

F (a) Move forward a step of length a > 0. The position of theturtle changes to (x′, y′, z′), where x′ = x + a�Hx, y′ = y + a�Hy

and z′ = z + a�Hz. A line segment is drawn between points(x, y, z) and (x′, y′, z′).

f(a) Move forward a step of length a without drawing a line.

+(a) Rotate around �U by an angle of a degrees. If a is positive, theturtle is turned to the left and if a is negative, the turn is tothe right.

&(a) Rotate around �L by an angle of a degrees. If a is positive,the turtle is pitched down and if a is negative, the turtle ispitched up.

/(a) Rotate around �H by an angle of a degrees. If a is positive, theturtle is rolled to the right and if a is negative, it is rolled tothe left.

It should be noted that symbols +, &, ∧, and / are used both asletters of the alphabet V and as operators in logical and arithmeticexpressions. Their meaning depends on the context.

The following examples illustrate the operation of parametric L-Row of treessystems with turtle interpretation. The first L-system is a coding ofa Koch construction generating a variant of the snowflake curve (Fig-ure 1.1 on page 2). The initiator (production predecessor) is the hy-potenuse AB of a right triangle ABC (Figure 1.36). The first and thefourth edge of the generator subdivide AB into segments AD and DB,while the remaining two edges traverse the altitude CD in opposite di-rections. From elementary geometry it follows that the lengths of thesesegments satisfy the equations

q = c − p and h =√

pq.

The edges of the generator can be associated with four triangles thatare similar to ABC and tile it without gaps. According to the relation-ship between curve construction by edge rewriting and planar tilings

1.10. Parametric L-systems 47

Figure 1.36: Construction of the generator for a “row of trees.” The edgesare associated with triangles indicated by ticks.

(Section 1.4.1), the generated curve will approximately fill the triangleABC. The corresponding L-system is given below:

#define c 1#define p 0.3#define q c − p#define h (p ∗ q) ∧ 0.5

ω : F (1)p1 : F (x) → F (x ∗ p) + F (x ∗ h) −−F (x ∗ h) + F (x ∗ q)

The resulting curve is shown in Figure 1.37a. In order to bettervisualize its structure, the angle increment has been set to 86◦ insteadof 90◦. The curve fills different areas with unequal density. This resultsfrom the fact that all edges, whether long or short, are replaced bythe generator in every derivation step. A modified curve that fills theunderlying triangle in a more uniform way is shown in Figure 1.37b. Itwas obtained by delaying the rewriting of shorter segments with respectto the longer ones, as specified by the following L-system.

ω : F (1, 0)p1 : F (x, t) : t = 0 → F (x ∗ p, 2) + F (x ∗ h, 1)−

−F (x ∗ h, 1) + F (x ∗ q, 0)p2 : F (x, t) : t > 0 → F (x, t − 1)

The next example makes use of node rewriting (Section 1.4.2). The Branchingstructureconstruction recursively subdivides a rectangular tile ABCD into two

tiles, AEFD and BCFE, similar to ABCD (Figure 1.38). The lengthsof the edges form the proportion

a

b=

b12a,

48 Chapter 1. Graphical modeling using L-systems

Figure 1.37: Two curves suggesting a “row of trees.” Curve (b) is from [95,page 57].

which implies that b = a/√

2. Each tile is associated with a single-pointframe lying in the tile center. The tiles are connected by a branchingline specified by the following L-system:

#define R 1.456ω : A(1)p1 : A(s) → F (s)[+A(s/R)][−A(s/R)]

(1.9)

The ratio of branch sizes R slightly exceeds the theoretical value of√2. As a result, the branching structure shown in Figure 1.39 is self-

avoiding. The angle increment was set arbitrarily to δ = 85◦.The L-system in equation (1.9) operates by appending segments

of decreasing length to the structures obtained in previous derivationsteps. Once a segment has been incorporated, its length does notchange. A structure with identical proportions can be obtained by

1.10. Parametric L-systems 49

Figure 1.38: Tiling associated with a space-filling branching pattern

Figure 1.39: A branching pattern generated by the L-system specified inequation (1.9)

50 Chapter 1. Graphical modeling using L-systems

Figure 1.40: Initial sequences of figures generated by the L-systems specifiedin equations (1.9) and (1.10)

appending segments of constant length and increasing the lengths ofpreviously created segments by constant R in each derivation step. Thecorresponding L-system is given below.

ω : Ap1 : A → F (1)[+A][−A]p2 : F (s) → F (s ∗ R)

(1.10)

The initial sequence of structures obtained by both L-systems arecompared in Figure 1.40. Sequence (a) emphasizes the fractal characterof the resulting structure. Sequence (b) suggests the growth of a tree.The next two chapters show that this is not a mere coincidence, andthe L-system specified in equation (1.10) is a simple, but in principlecorrect, developmental model of a sympodial branching pattern foundin many herbaceous plants and trees.

Chapter 2

Modeling of trees

Computer simulation of branching patterns has a relatively long history. Cellular–spacemodelsThe first model was proposed by Ulam [149], (see also [138, pages 127–

131]), and employed the concept of cellular automata that had beendeveloped earlier by von Neumann and Ulam [156]. The branchingpattern is created iteratively, starting with a single colored cell in atriangular grid, then coloring cells that touch one and only one vertexof a cell colored in the previous iteration step.

This basic idea gave rise to several extensions. Meinhardt [97, Chap-ter 15] substituted the triangular grid with a square one, and used theresulting cellular space to examine biological hypotheses related to theformation of net-like structures. In addition to pure branching patterns,his models capture the effect of branch reconnection or anastomosis thatmay take place between the veins of a leaf. Greene [54] extended cellu-lar automata to three dimensions, and applied the resulting voxel spaceautomata to simulate growth processes that react to the environment.For instance, Figure 2.1 presents the growth of a vine over a house.Cohen [15] simulated the development of a branching pattern using ex-pansion rules that operate in a continuous “density field” rather thana discrete cellular or voxel space.

The common feature of these approaches is the emphasis on inter-actions between various elements of a growing structure, as well as thestructure and the environment. Although interactions clearly influencethe development of real plants, they also add to the complexity of themodels. This may explain why simpler models, ignoring even such fun-damental factors as collisions between branches, have been prevalent todate. The first model in that category was proposed by Honda [65] who Honda’s modelstudied the form of trees using the following assumptions (Figure 2.2).

• Tree segments are straight and their girth is not considered.

• A mother segment produces two daughter segments through onebranching process.

52 Modeling of trees

Figure 2.1: Organic architecture by Greene [54]

• The lengths of the two daughter segments are shortened by con-stant ratios, r1 and r2, with respect to the mother segment.

• The mother segment and its two daughter segments are containedin the same branch plane. The daughter segments form constantbranching angles, a1 and a2, with respect to the mother branch.

• The branch plane is fixed with respect to the direction of gravityso as to be closest to a horizontal plane.1 An exception is madefor branches attached to the main trunk. In this case, a constantdivergence angle α between consecutively issued lateral segmentsis maintained (cf. Chapter 4).

By changing numerical parameters, Honda obtained a wide vari-ety of tree-like shapes. With some improvements [38], his model wasapplied to investigate branching patterns of real trees [39, 66, 67, 68].Subsequently, different rules for branching angles were proposed to cap-ture the structure of trees in which planes of successive bifurcations areperpendicular to each other [69]. The results of Honda served as a basisfor the tree models proposed by Aono and Kunii [2]. They suggested

1More formally, the line perpendicular to the mother segment and lying in thebranch plane is horizontal.

Modeling of trees 53

Figure 2.2: Specification of tree geometry according to Honda [65]

several extensions to his model, the most important of which was thebiasing of branch positions in a particular direction, applied to producethe effects of wind, phototropism and gravity. A similar concept wasintroduced previously by Cohen [15], while more accurate physically-based methods for branch bending were developed by de Reffye [28]and Armstrong [4].

The models of Honda and Aono and Kunii were rendered using Realismstraight lines of constant or varying width to represent “tree skele-tons.” A substantial improvement in the realism of synthetic imageswas achieved by Bloomenthal [11] and Oppenheimer [105], who intro-duced curved branches, carefully modeled surfaces around branchingpoints, and textured bark and leaves (Figure 2.3).

The approaches stemming from the work of Honda defined branching Stochasticmodelsstructures using deterministic algorithms. In contrast, stochastic mech-

anisms are essential to the group of tree models proposed by Reeves andBlau [119], de Reffye et al. [30], and Remphrey, Neal and Steeves [120].Although these models were described using different terminologies,they share the basic paradigm of specifying tree structures in terms ofprobabilities with which branches are formed. The work of Reeves andBlau aimed at producing tree-like shapes without delving into biologicaldetails of the modeled structures (Figure 2.4). In contrast, de Reffyeet al. [29] used a stochastic approach to simulate the development of Approach of

de Reffyereal plants by modeling the activity of buds at discrete time intervals.Given a clock signal, a bud can either:

• do nothing,

• become a flower,

• become an internode terminated by a new straight apex and oneor more lateral apices subtended by leaves, or

• die and disappear.

54 Modeling of trees

Figure 2.3: Acer graphics by Bloomenthal [11]

Figure 2.4: A forest scene by Reeves [119] c©1984 Pixar

Modeling of trees 55

Figure 2.5: Oil palm tree canopy from CIRAD Modelisation Laboratory

These events occur according to stochastic laws characteristic for eachvariety and each species. The geometric parameters, such as the lengthand diameter of an internode, as well as branching angles, are alsocalculated according to stochastic laws.

The basic types of developmental rules incorporated into this methodcorrespond to the 23 different types of tree architectures identified byHalle, Oldeman and Tomlinson [58]. Detailed models of selected plantspecies were developed and are described in the literature [16, 20, 26,27, 76]. A sample tree model is shown in Figure 2.5. The approach ofRemphrey [120, 121, 122] is similar to that of de Reffye, except thatlarger time steps are used (one year in the model of bearberry describedin [120]). Consequently, the stochastic rules must describe the entireconfiguration of lateral shoots that can be formed over a one-year pe-riod.

The application of L-systems to the generation of botanical trees was Application ofL-systemsfirst considered by Aono and Kunii [2]. They referred to the original

definition of L-systems [82] and found them unsuitable to model thecomplex branching patterns of higher plants. However, their argumentsdo not extend to parametric L-systems with turtle interpretation. Forexample, the L-system in Figure 2.6 implements those tree models ofHonda [65] in which one of the branching angles is equal to 0, yieldinga monopodial structure with clearly delineated main and lateral axes(see Section 3.2 for a formal characterization).

56 Modeling of trees

n = 10#define r1 0.9 /* contraction ratio for the trunk */#define r2 0.6 /* contraction ratio for branches */#define a0 45 /* branching angle from the trunk */#define a2 45 /* branching angle for lateral axes */#define d 137.5 /* divergence angle */#define wr 0.707 /* width decrease rate */

ω : A(1,10)p1: A(l,w) : *→ !(w)F(l)[&(a0)B(l*r2,w*wr)]/(d)A(l*r1,w*wr)p2: B(l,w) : *→ !(w)F(l)[-(a2)$C(l*r2,w*wr)]C(l*r1,w*wr)p3: C(l,w) : *→ !(w)F(l)[+(a2)$B(l*r2,w*wr)]B(l*r1,w*wr)

Figure 2.6: Examples of the monopodial tree-like structures of Honda [65],generated using L-systems

Modeling of trees 57

According to production p1, the apex of the main axis A produces Monopodialbranchingan internode F and a lateral apex B in each derivation step. Con-

stants r1 and r2 specify contraction ratios for the straight and lateralsegments, a0 and a2 are branching angles and d is the divergence angle.Module !(w) sets the line width to w, thus production p1 decreases Branch widththe width of the daughter segments with respect to the mother seg-ment by the factor wr = 0.707. This constant satisfies a postulate byLeonardo da Vinci [95, page 156], according to which “all the branchesof a tree at every stage of its height when put together are equal inthickness to the trunk below them.” In the case where a mother branchof diameter w1 gives rise to two daughter branches of equal diameterw2, this postulate yields the equation w2

1 = 2w22, which gives a value for

wr equal to w2/w1 = 1/√

2 ≈ 0.707. A general discussion of the rela-tionships between the diameters of the mother and daughter branchesis included in the book by Macdonald [94, pages 131–135].

Productions p2 and p3 describe subsequent development of the lat-eral branches. In each derivation step, the straight apex (either B or C)issues a lateral apex of the next order at angle a2 or −a2 with respectto the mother axis. Two productions are employed to create lateralapices alternately to the left and right. The symbol $ rolls the turtle Keeping

turtle’sorientation

around its own axis so that vector �L pointing to the left of the turtle(Section 1.5) is brought to a horizontal position. Consequently, thebranch plane is “closest to a horizontal plane,” as required by Honda’smodel. From a technical point of view, $ modifies the turtle orientationin space according to the formulae

�L =�V × �H

|�V × �H|and �U = �H × �L,

where vectors �H, �L and �U are the heading, left and up vectors associatedwith the turtle, �V is the direction opposite to gravity, and |�A| denotesthe length of vector �A. The tree-like structures shown in Figure 2.6were generated using the values of constants listed in Table 2.1, andcoincide with the structures presented by Honda.

Figure r1 r2 a0 a2

a 0.9 0.6 45 45b 0.9 0.9 45 45c 0.9 0.8 45 45d 0.9 0.7 30 -30

Table 2.1: Constants for the monopodial tree structures in Figure 2.6

58 Modeling of trees

A slightly different L-system, specified in Figure 2.7, captures sym-Sympodialbranching podial structures, where both daughter segments form non-zero angles

with the mother segment. In this case the activity of the main apex isreduced to the formation of the trunk F and a pair of lateral apices B(production p1). The subsequent branching pattern is captured by pro-duction p2. The sample structures in Figure 2.7 were obtained using theconstants listed in Table 2.2, and correspond to the models presentedby Aono and Kunii.

The L-systems considered so far have been designed in a mannerthat emphasizes their relation to the models described in the litera-ture. Specifically, all segments are assigned their final length at thetime of creation, and no further elongation occurs. As pointed out inSection 1.10.3, similar structures can be obtained by creating new seg-ments of constant length and increasing the lengths of previously cre-ated segments by a constant factor in each derivation step. A sampleL-system constructed according to this paradigm is given in Figure 2.8.

The overall structure of the tree is defined by production p1. InTernarybranching each derivation step, apex A produces three new branches terminated

by their own apices. Parameter w and constant vr relate the width ofthe mother branch w1 to that of the daughter branches w2. Accordingto da Vinci’s postulate w2

1 = 3w22, thus vr = w1/w2 =

√3 ≈ 1.732.

Productions p2 and p3 capture the gradual elongation of branches andthe increase in their diameter over time.

The bending of branches is simulated by slightly rotating the tur-Tropismtle in the direction of a predefined tropism vector �T after drawing eachsegment (Figure 2.9). The orientation adjustment α is calculated us-ing the formula α = e |�H × �T |, where e is a parameter capturing axissusceptibility to bending. This heuristic formula has a physical moti-vation; if �T is interpreted as a force applied to the endpoint of vector �H,and �H can rotate around its starting point, the torque is equal to �H ×�T .The parameters relevant to the generation of the tree-like structures inFigure 2.8 are listed in Table 2.3. A more realistic rendering of the treein Figure 2.8d is presented in Figure 2.10.

Figure r1 r2 a1 a2

a 0.9 0.7 5 65b 0.9 0.7 10 60c 0.9 0.8 20 50d 0.9 0.8 35 35

Table 2.2: Constants for the sympodial tree structures in Figure 2.7

Modeling of trees 59

n = 10#define r1 0.9 /* contraction ratio 1 */#define r2 0.7 /* contraction ratio 2 */#define a1 10 /* branching angle 1 */#define a2 60 /* branching angle 2 */#define wr 0.707 /* width decrease rate */

ω : A(1,10)p1 : A(l,w) : *→ !(w)F(l)[&(a1)B(l*r1,w*wr)]

/(180)[&(a2)B(l*r2,w*wr)]p2 : B(l,w) : *→ !(w)F(l)[+(a1)$B(l*r1,w*wr)]

[-(a2)$B(l*r2,w*wr)]

Figure 2.7: Examples of the sympodial tree-like structures of Aono andKunii [2], generated using L-systems

60 Modeling of trees

#define d1 94.74 /* divergence angle 1 */#define d2 132.63 /* divergence angle 2 */#define a 18.95 /* branching angle */#define lr 1.109 /* elongation rate */#define vr 1.732 /* width increase rate */

ω : !(1)F(200)/(45)Ap1 : A : * → !(vr)F(50)[&(a)F(50)A]/(d1)

[&(a)F(50)A]/(d2)[&(a)F(50)A]p2 : F(l) : *→ F(l*lr)p3 : !(w) : *→ !(w*vr)

Figure 2.8: Examples of tree-like structures with ternary branching

Modeling of trees 61

Figure 2.9: Correction α of segment orientation �H due to tropism �T

Figure d1 d2 a lr �T e na 94.74 132.63 18.95 1.109 0.00,-1.00,0.00 0.22 6b 137.50 137.50 18.95 1.109 0.00,-1.00,0.00 0.14 8c 112.50 157.50 22.50 1.790 -0.02,-1.00,0.00 0.27 8d 180.00 252.00 36.00 1.070 -0.61,0.77,-0.19 0.40 6

Table 2.3: Constants for the tree structures in Figure 2.8

Figure 2.10: Medicine Lake by Musgrave et al. [101]

62 Modeling of trees

Figure 2.11: A surrealistic elevator

The examples given above demonstrate that the tree models of Honda,Conclusionsas well as their derivatives studied by Aono and Kunii, can be expressedusing the formalism of L-systems. In a separate study, Shebell [130] alsoshowed that L-systems can be applied to generate the architectural treemodels of Halle, Oldeman and Tomlinson [58]. These results indicatethat L-systems may play an important role as a tool for biologically-correct simulation of tree development and synthesis of realistic treeimages. However, the tree-like shapes created so far are rather generic(Figure 2.11), and models of particular tree species, directly based onbiological data, are yet to be developed. L-systems have found moreapplications in the domain of realistic modeling of herbaceous plants,discussed in the next chapter.

Chapter 3

Developmental models ofherbaceous plants

The examples of trees presented in the previous chapter introduce L-systems as a plant modeling tool. They also illustrate one of the moststriking features of the generative approach to modeling, called database amplification [136]. This term refers to the generation of complex-looking objects from very concise descriptions – in our case, L-systemscomprised of small numbers of productions. Yet in spite of the smallsize, the specification of L-systems is not a trivial task.

In the case of highly self-similar structures, the synthesis methodsbased on edge rewriting and node rewriting are of assistance, as illus-trated by the examples considered in Section 1.10.3. However, a moregeneral approach is needed to model the large variety of developmentalpatterns and structures found in nature.

The methodology presented in this chapter is based on the simu- Developmentalmodelslation of the development of real plants. Thus, a particular form is

modeled by capturing the essence of the developmental process thatleads to this form. This approach has two distinctive features.

• Emphasis on the space-time relation between plant parts.In many plants, organs in various stages of development can beobserved at the same time. For example, some flowers may stillbe in the bud stage, others may be fully developed, and stillothers may have been transformed into fruits. If developmentis simulated down to the level of individual organs, such phaseeffects are reproduced in a natural way.

• Inherent capability of growth simulation. The mathemati-cal model can be used to generate biologically correct images ofplants at different ages and to create sequences of images illus-trating plant development in time.

64 Chapter 3. Developmental models

The models are constructed under the assumption that organismscontrol the important aspects of their own development. Accordingto Apter [3, page 44], this simplification must be accepted as a nec-essary evil, as long as the scope of the mathematical model is limitedto an isolated plant. Consequently, this chapter focuses on the mod-eling and generation of growth sequences of herbaceous or non-woodyHerbaceous

plants plants, since internal control mechanisms play a predominant role intheir development. In contrast, the form of woody plants is determinedto a large extent by the environment, competition among branches andtrees, and accidents [164].

3.1 Levels of model specification

L-systems can be constructed with a variety of objectives in mind, rang-ing from a general classification of branching structures to detailed mod-els suitable for image synthesis purposes. Accordingly, the L-systemspresented in this chapter are specified at three levels of detail. The mostabstract level, called partial L-systems, employs the notation of nonde-Partial

L-systems terministic OL-systems to define the realm of possibilities within whichstructures of a given type may develop. Partial L-systems capture themain traits characterizing structural types, and provide a formal basisfor their classification. Control mechanisms that resolve nondetermin-ism are introduced in the next level, termed L-system schemata.1 TheL-system

schemata topology of individual plants and temporal aspects of their develop-ment are described at this level. Schemata are of particular interestfrom a biological point of view, as they provide an insight into themechanisms that control plant development in nature. The geomet-ric aspects are added in complete L-systems that include informationComplete

L-systems concerning growth rates of internodes, the values of branching angles,and the appearance of organs. The difference between all three levels isillustrated using models of a single-flower shoot as a running example.

3.1.1 Partial L-systems

Consider the development of a shoot which, after a period of vegetativeSingle-flowershoot growth, produces a single flower. The partial L-system is given below.

ω : ap1 : a → I[L]ap2 : a → I[L]Ap3 : A → K

(3.1)

The lower-case symbol a represents the vegetative apex, while the upper-case A is the flowering apex, capable of forming reproductive organs.

1In the literature, the term “scheme” is also used to denote the class of L-systemswith the same alphabet and productions, but with different axioms [62, page 54].

3.1. Levels of model specification 65

Figure 3.1: Single-flower shoot

A derivation step corresponds to a plastochron, defined as the timeinterval between the production of successive internodes by the apex.At each step apex a has a choice of forming either leaf L, internode Iand new apex a (production p1), or forming the same structures andtransforming itself into a flowering apex A (p2), which subsequentlycreates flower K (p3). Once this transformation or developmental switchhas taken place it cannot be reversed, since there is no rule allowingthe transformation of A to a. Examples of strings generated by theL-system specified in equation (3.1) are given below.

a a aI[L]A I[L]a I[L]aI[L]K I[L]I[L]A I[L]I[L]aI[L]K I[L]I[L]K I[L]I[L]I[L]A

I[L]I[L]K I[L]I[L]I[L]KI[L]I[L]I[L]K

A diagrammatic representation of a single-flower inflorescence is shownin Figure 3.1.

3.1.2 Control mechanisms in plants

A partial L-system does not specify the moments in which develop- L-systemschematamental switches occur. The timing of these switches is specified at the

level of L-system schemata, which incorporate mechanisms that con-trol plant development. In biology, these mechanisms are divided intotwo classes depending on the way information is transferred between Lineage vs.

interactionmodules. The term lineage (or cellular descent) refers to the transferof information from an ancestor cell or module to its descendants. Incontrast, interaction is the mechanism of information exchange betweenneighboring cells (for example, in the form of nutrients or hormones).Within the formalism of L-systems, lineage mechanisms are represented

66 Chapter 3. Developmental models

by context-free productions found in OL-systems, while the simulationof interaction requires the use of context-sensitive 1L-systems and 2L-systems.2 Several specific mechanisms are listed below. Although theyare described from the modeling perspective, a relation to physiologicalprocesses observed in nature can often be found.

Stochastic mechanism

The simplest method for implementing a developmental switch is touse a stochastic L-system. In this case the vegetative apex a has aprobability π1 of staying in the vegetative state, and π2 of transformingitself into a flowering apex A.

ω : ap1 : a

π1→ I[L]ap2 : a

π2→ I[L]A

p3 : A1→ K

The probability distribution (π1, π2) is found experimentally, with π1 +π2 = 1.

The effect of environment

Many plants change from a vegetative to a flowering state in responseto environmental factors such as temperature or the number of day-light hours. Such effects can be modeled using one set of productions(called a table) for some number of derivation steps, then replacing itTable

L-systems by another set:

Table 1 Table 2ω : a p1 : a → I[L]Ap1 : a → I[L]a p2 : A → K

The concept of table L-systems (TOL-systems) was introduced andformalized by Rozenberg [62, 127]. Note that the use of tables providesonly a partial solution to the problem of specifying the switching time,since a control mechanism external to the L-system is needed to selectthe appropriate table.

Delay mechanism

The delay mechanism operates under the assumption that the apexundergoes a series of state changes that postpone the switch until aparticular state is reached.

2The clarity of this dichotomy is somewhat obscured by parametric OL-systems,which can simulate the operation of context-sensitive L-systems using an infiniteset of parameter values.

3.1. Levels of model specification 67

This is captured by the following L-system in the case of a single-flowershoot.

ω : a0

pi : ai → I[L]ai+1 0 ≤ i ≤ n − 1pn : an → I[L]Apn+1 : A → K

According to this model, the apex counts the leaves it produces. Whileit may seem strange that a plant would count, it is known that someplant species produce a fixed number of leaves before they start flow-ering.

Accumulation of components

A developmental mechanism based on the accumulation of componentsis similar to that of delay, but emphasizes the physiological nature ofthe counting process. According to this approach, counting is achievedby a monotonic increase or decrease in the concentration of certain cellcomponents. This process can be captured by the following parametricL-system:

ω : a(0)p1 : a(c) : c < C → I[L]a(c + ∆c)p2 : a(c) : c ≥ C → I[L]Ap3 : A : ∗ → K

(3.2)

The parameter c indicates current concentration of the controlling com-ponents in the apex a. In each derivation step, this concentration isincreased by a constant ∆c. The developmental switch occurs whenthe concentration reaches the threshold value C.

Development controlled by a signal

In many plants, the switch from a vegetative to a flowering state iscaused by a flower-inducing signal transported from the basal leavestowards the apex. The time of signal initiation is determined usingone of the previously described methods, for example by counting. Asample L-system is given below.

ω : D(1)a(1)p1 : a(i) : i < m → a(i + 1)p2 : a(i) : i = m → I[L]a(1)p3 : D(i) : i < d → D(i + 1)p4 : D(i) : i = d → S(1)p5 : S(i) : i < u → S(i + 1)p6 : S(i) : i = u → εp7 : S(i) < I : i = u → IS(1)p8 : S(i) < a(j) : ∗ → I[L]Ap9 : A : ∗ → K

68 Chapter 3. Developmental models

The apex a produces internodes I and leaves L on the main axis (p2).The time between the production of two consecutive segments, or theplastochron of the main axis, is equal to m derivation steps (p1). Aftera delay of d steps (p3), a signal S is sent from the plant base towardsthe apices (p4). This signal is transported along the main axis witha delay of u steps per internode I (p5,p7). Production p6 removes thesignal from a node after it has been transported along the structure(ε stands for the empty string). When the signal reaches the apex, a istransformed into flowering state A (p8) which yields flower K (p9). Notethat the signal has to propagate faster than one node per plastochron(u < m), otherwise it would not be able to catch up with the apex. Theabove processes are illustrated by the following developmental sequence,for d = 4, m = 2 and u = 1.

D(1)a(1)D(2)a(2)D(3)I[L]a(1)D(4)I[L]a(2)S(1)I[L]I[L]a(1)IS(1)[L]I[L]a(2)I[L]IS(1)[L]I[L]a(1)I[L]I[L]IS(1)[L]a(2)I[L]I[L]I[L]AI[L]I[L]I[L]K

Although the above model may appear unnecessarily complicated, sig-nals are indispensable in the simulation of complex flowering sequencesdiscussed later.

3.1.3 Complete models

The L-systems considered so far are not directly suitable for image syn-thesis purposes. To this end, they must be completed with geometricinformation. The relation between an L-system scheme and a corre-sponding complete L-system is discussed using the model of crocusesshown in Figure 3.2 as an example.

The development is controlled using a delay expressed as an accu-Crocusmulation mechanism (equation (3.2)). In contrast to L-system schemesin which symbols represent module types, the L-system in Figure 3.2is specified in terms of turtle symbols. Production p1 describes the cre-ation of successive internodes F and leaves L by the vegetative apex a.The leaves branch from the stem at an angle of 30◦ and spiral aroundthe main axis with a divergence angle equal to 137.5◦ (see Chapter 4).Productions p2 and p3 describe the developmental switch and the cre-ation of flower K taking place respectively in steps Ta and Ta+1. Pro-ductions p4 and p5 capture the development of leaves and flowers untilthey reach their final shapes TL and TK steps after creation. For each

3.1. Levels of model specification 69

#define Ta 7 /* developmental switch time */#define TL 9 /* leaf growth limit */#define TK 5 /* flower growth limit */#include L(0),L(1),...,L(TL) /* leaf shapes */#include K(0),K(1),...,K(TK) /* flower shapes */

ω : a(1)p1 : a(t) : t<Ta → F(1)[&(30)∼L(0)]/(137.5)a(t+1)p2 : a(t) : t=Ta → F(20)Ap3 : A : * → ∼K(0)p4 : L(t) : t<TL → L(t+1)p5 : K(t) : t<TK → K(t+1)p6 : F(l) : l<2 → F(l+0.2)

Figure 3.2: Crocuses

70 Chapter 3. Developmental models

branchmain apex main apexterminates continues

all some all somelateral lateral lateral lateralapices apices apices apices

terminate continue terminate continueterminal sympodial monopodial polypodial

Table 3.1: Basic growth patterns of branching structures

value of parameter t, the corresponding organ shapes L(t) and K(t) aremodeled using bicubic patches incorporated into the plant structure asdescribed in Section 5.1. Production p6 specifies the gradual elongationof internodes.

3.2 Branching patterns

The partial L-system in equation (3.1) and the related schemata employSubapicalgrowth subapical growth mechanisms in which new branches are created exclu-

sively by apices. All herbaceous plants develop this way. The archi-tecture of a branching structure is to a large extent determined by therelationships between terminal and continuing apices. While a contin-uing apex produces branches again and again, a terminal one eithergives rise to an appendage such as a flower or dies. The possible com-binations are listed in Table 3.1. Of the four terms assigned to thesepossibilities, two are commonly used in biology, namely, sympodial andmonopodial, while the other two terms are introduced here to denotethe remaining cases, usually not characterized in the literature.

The above branching patterns can be represented conveniently usingExpressionusingL-systems

partial bracketed L-systems. Let A,B,C denote continuing apices, Xa terminal apex, and I an internode. The terminal and sympodialpatterns are characterized by rules of the form

A → I[B]n[X]mX,

with n = 0, m > 0 in the case of terminal patterns and n > 0, m ≥ 0in the case of sympodial patterns. The important property is that themain apex terminates its development in all these cases.

Monopodial and polypodial patterns have rules of the form

A → I[B]n[X]mC,

with n = 0, m > 0 in the case of monopodial patterns and n > 0,m ≥ 0 in the case of polypodial ones. In these cases, it is importantthat the apex remains active.

3.3. Models of inflorescences 71

The terms defined here apply to branching structures in general,whether they only produce vegetative organs (branches and leaves) orreproductive organs as well. In the latter case, the terminal organsdevelop from flower buds to flowers to fruits but do not give rise tovegetative structures (there are exceptions to this statement but theycan be neglected here).

3.3 Models of inflorescences

The following discussion focuses on the modeling of compound floweringstructures or inflorescences. In some cases an entire shoot system canbe considered an inflorescence, in others only some of the branchesbear flowers and are inflorescences. Inflorescence architecture is anelaboration of branching structures in general.

In the domain of botanical applications of L-systems, the study Classificationof inflorescences has played a particularly visible role [44, 45, 46, 47,77, 86]. Unfortunately, the terms used for the various inflorescencetypes are not uniform in the literature. Besides a purely morphologicalterminology, attempts have been made to construct a “typological”terminology, expressing the “essential” features of flowering structures[144, 145, 157, 158]. However, these terms are not generally accepted[12]. A compromise has been proposed by D. and U. Muller-Doblies[100], which serves as a basis for the classification that guides thispresentation.

3.3.1 Monopodial inflorescences

Simple racemes (open)

Racemes are characteristically monopodial inflorescences; a shoot haslateral apices with terminal structures and a main apex that continuesto development. A raceme is open if the main apex does not form aflower. The partial L-system for this widely occurring type of inflores-cence is:

ω : ap1 : a → I[L]ap2 : a → I[L]Ap3 : A → I[K]A

(3.3)

This system differs from that modeling a shoot with a single flower(equation (3.1)) only in production p3. Here it is designed to repeatedlyproduce lateral flowers (Figure 3.3), while in the previous system Aproduces a single flower.

72 Chapter 3. Developmental models

a b

Figure 3.3: Open racemes: (a) elongated form, (b) planar form

Figure 3.4: Lily-of-the-valley

3.3. Models of inflorescences 73

The flowering sequence in open racemes is always acropetal (from Acropetalsequencebase to top). This can be observed after substituting production p3 in

the L-system in equation (3.3) with productions p′3 and p4, which use

indexed symbols Ki to denote subsequent stages of flower development.

p′3 : A → I[IK0]A

p4 : Ki → Ki+1, i ≥ 0

The indexed notation Ki → Ki+1 stands for a (potentially infinite)set of productions K0 → K1, K1 → K2, K2 → K3,.... The developmen-tal sequence begins as follows:

AI[IK0]AI[IK1]I[IK0]AI[IK2]I[IK1]I[IK0]AI[IK3]I[IK2]I[IK1]I[IK0]A· · ·

At each developmental stage the inflorescence contains a sequence Lily-of-the-valleyof flowers of different ages. The flowers newly created by the apex are

delayed in their development with respect to the older ones situated atthe stem base. Graphically, this effect is illustrated by the model ofa lily-of-the-valley shown in Figure 3.4. The following quotation fromd’Arcy Thompson [143] applies:

A flowering spray of lily-of-the-valley exemplifies a growth-gradient, after a simple fashion of its own. Along the stalkthe growth-rate falls away; the florets are of descending age,from flower to bud; their graded differences of age lead toan exquisite gradation of size and form; the time-intervalbetween one and another, or the “space-time relation” be-tween them all, gives a peculiar quality – we may call itphase-beauty – to the whole.

Another example of “phase beauty” can be seen in the shoot of shep- Capsellaherd’s purse (Capsella bursa-pastoris) shown in Figure 3.5. Productionsp1, p2 and p3 describe the activities of the apex in the vegetative andflowering states, in accordance with the L-system in equation (3.3). Thedevelopmental switch is implemented using a delay mechanism. Pro-ductions p4 and p5 capture the linear elongation of internodes in time,while p6 and p7 describe the gradual increase of the angle at whichthe flower stalks branch from the main stem. Productions p8, p9 andp11 specify the shapes of leaves L, flower petals K and fruits X usingdevelopmental surface models discussed in Section 5.2. Production p10

controls the flowering time. Symbol % in the successor of productionp11 simulates the fall of petals by cutting them off the structure atthe time of fruit formation. The default value of the angle incrementcorresponding to the symbol + with no parameter is 18◦.

74 Chapter 3. Developmental models

ω : I(9)a(13)p1 : a(t) : t>0 → [&(70)L]/(137.5)I(10)a(t-1)p2 : a(t) : t=0 → [&(70)L]/(137.5)I(10)Ap3 : A : * → [&(18)u(4)FFI(10)I(5)X(5)KKKK]

/(137.5)I(8)Ap4 : I(t) : t>0 → FI(t-1)p5 : I(t) : t=0 → Fp6 : u(t) : t>0 → &(9)u(t-1)p7 : u(t) : t=0 → &(9)p8 : L : * → [{.-FI(7)+FI(7)+FI(7)}]

[{.+FI(7)-FI(7)-FI(7)}]p9 : K : * → [&{.+FI(2)--FI(2)}]

[&{.-FI(2)++FI(2)}]/(90)p10 : X(t) : t>0 → X(t-1)p11 : X(t) : t=0 → ∧(50)[[-GGGG++[GGG[++G{.].].].

++GGGG.--GGG.--G.}]%

Figure 3.5: Development of Capsella bursa-pastoris. Every fourth derivationstep is shown.

3.3. Models of inflorescences 75

Figure 3.6: Apple twig

Simple raceme (closed)

The inflorescence of an apple tree (Figure 3.6) provides an example ofa closed raceme. In this case, the main apex eventually terminates itsdevelopment and produces a terminal flower (Figure 3.7). The corre-sponding partial L-system is given below.

ω : ap1 : a → I[L]ap2 : a → I[L]Ap3 : A → I[K]Ap4 : A → K

Developmental switches are associated with two symbols, a and A.Thus, in order to obtain an L-system scheme it is necessary to specifyhow both of these switches will be controlled.

The flowering sequence is usually acropetal but could also be basi-petal, i.e., progressing downward after the formation of the terminalflower on the main axis. In the latter case a basipetal signal, as dis-cussed in Section 1.8, can be applied to induce the transformation ofdormant flower buds into flowers.

76 Chapter 3. Developmental models

a b

Figure 3.7: Closed racemes: (a) elongated form, (b) planar form

Compound raceme (open dibotryoid)

Racemes can also occur on complex branching structures. The sim-plest of these inflorescences is one with open racemes on the first orderbranches as well as on the main axis (Figure 3.8a). This two-level com-pound structure (thus dibotryoid) is described by the following partialL-system.

ω : ap1 : a → I[L]ap2 : a → I[L]Ap3 : A → I[L][b]Ap4 : A → I[L][b]Bp5 : b → I[L]bp6 : b → I[L]Bp7 : B → I[K]B

(3.4)

Three developmental transformations are necessary: the first for thechange from leaf to branch creation along the main axis (production p2),the second for the switch from branching to lateral flower creation onthe main axis (p4), and the third for the transition from leaf to lateralflower formation along the first-order branches (p6). Each branch issubtended by a leaf, which is why productions p3 and p4 specify twoappendages L and b. Branches with flowers K need not have subtendingleaves, which is reflected in production p7.

Within each component raceme, the flowering sequence is alwaysSingle-signalmodel acropetal, but the timing of switches has a crucial impact on the over-

all flowering sequence and appearance of the plant. For example, let usassume that the switch from leaf to branch production is controlled bya delay, while the remaining two switches are caused by an acropetalflower-inducing signal (representing the hormone florigen). Such a de-velopment is captured by L-system 3.1 (see below). Initially, the veg-

3.3. Models of inflorescences 77

a b

Figure 3.8: Dibotryoids: (a) open, (b) closed

#define d 13 /* delay for sending florigen */#define e 3 /* delay for creating branches */#define m 2 /* plastochron - main axis */#define n 3 /* plastochron - lateral axis */#define u 1 /* signal delay - main axis */#define v 1 /* signal delay - lateral axis */

ω : S(0)a(1,0)p1 : a(t,c) : t<m → a(t+1,c)p2 : a(t,c) : (t=m)&(c<e)→ I(0,u)[L]a(1,c+1)p3 : a(t,c) : (t=m)&(c=e)→ I(0,u)[L][b(1)]a(1,c)p4 : b(t) : t<n → b(t+1)p5 : b(t) : t=n → I(0,v)[L]b(1)p6 : S(t) : * → S(t+1)p7 : S(t) < I(i,j) : t=d → I(1,j)p8 : I(i,j) : (0<i)&(i<j)→ I(i+1,j)p9 : I(i,j) < I(k,l) : (i=j)&(k=0)→ I(1,l)p10: I(i,j) < a(k,l) : i>0 → I[L][b(1)]Bp11: I(i,j) < b(k) : i>0 → I[L]Bp12: B : * → I[K]B

L-system 3.1: A model of dibotryoids

78 Chapter 3. Developmental models

etative apex a creates internodes I and leaves L with plastochron m(productions p1 and p2). After the creation of e leaves a developmen-tal switch occurs, and apex a starts creating branches with the sameplastochron (p3). The change of state is indicated by the value of thesecond parameter in the module a(t, c), which is now equal to e. Thelateral apices b create internodes and leaves with plastochron n (p4 andp5). After a delay of d steps from the beginning of the simulation (p6),the flowering signal is introduced to the basal internode (p7), as indi-cated by a non-zero value of the first parameter in the module I(i, j).The signal is passed along an axis at the rate of j steps per internode(p8 and p9), where j = u for the main axis and j = v for the lateralaxes. These rates are assigned to internodes by productions p2, p3 andp5. When the signal reaches an apex (either a or b), the apex is trans-formed into flowering state B (p10 and p11). From then on, new flowersK are produced in each derivation step (p12).

In order to analyze the plant structure and flowering sequence re-Model analysissulting from the above development, let Tk denote the time at whichapex b of the k-th lateral axis is transformed into the flowering state,and lk denote the length of this axis (expressed as the number of in-ternodes) at the transformation time. It is assumed here that the first eleaves count as lateral axes, thus k > e. Since it takes km time units toproduce k internodes along the main axis and lkn time units to producelk internodes on the lateral axis, we obtain:

Tk = km + lkn

On the other hand, the transformation occurs when the signal reachesthe apex. The signal is sent d time units after the development starts.It uses ku time units to travel through k zero-order internodes and lkvtime units to travel through lk first-order internodes:

Tk = d + ku + lkv

Solving the above system of equations for lk and Tk (and ignoring forsimplicity some inaccuracy due to the fact that this system does notguarantee integer solutions), we obtain:

Tk = kun − vm

n − v+ d

n

n − v

lk = −km − u

n − v+

d

n − v

In order to analyze the above solutions, let us first notice that the signaltransportation delay v must be less than the plastochron of the lateralaxes n, otherwise the signal would never reach the lateral apices. Underthis assumption, the sign of the expression ∆ = un−vm determines theoverall flowering sequence, which is acropetal for ∆ > 0 (Figure 3.9)and basipetal for ∆ < 0 (Figure 3.10). If ∆ = 0, all flowering switches

3.3. Models of inflorescences 79

Figure 3.9: An acropetal flowering sequence in an open dibotryoid: m = 2,n = 3, u = v = 1, ∆ = 0.5; derivation lengths: 15−18−21−24−27−30−33

occur simultaneously. The sign of the expression m − u determineswhether the vegetative part of the shoot is more developed at the base(m − u < 0) or near the top of the structure (m − u > 0). Figure 3.11shows a model of a member of the mint family that exhibits a basipetalflowering sequence.

Compound racemes (closed dibotryoids)

This inflorescence type differs from the previous one only in that eachbranch, including the main axis, bears a terminal flower (Figure 3.8b).A partial L-system can be obtained from that of equation (3.4) byadding one more production:

p8 : B → K

80 Chapter 3. Developmental models

Figure 3.10: A basipetal flowering sequence in an open dibotryoid: m = 2,n = 5, u = 1, v = 3, ∆ = −0.5; derivation lengths: 16− 20− 24− 28− 32−36 − 40

3.3. Models of inflorescences 81

Figure 3.11: A mint

Compound raceme (closed tribotryoid)

Racemic inflorescences can be compounded to a higher number of lev-els. The following is a partial L-system for a closed tribotryoid inflo-rescence, where closed racemes occur on second-order branches as wellas on the terminal portions of first-order branches and of the main axis(Figure 3.12). The developmental process involves six developmentaltransformations.

ω : ap1 : a → I[L]ap2 : a → I[L]Ap3 : A→ I[L][b]Ap4 : A→ I[L][b]Bp5 : b → I[L]bp6 : b → I[L]B

p7 : B→ I[L][c]Bp8 : B→ I[L][c]Cp9 : c → I[L]cp10 : c → I[L]Cp11 : C→ I[K]Cp12 : C→ K

82 Chapter 3. Developmental models

Figure 3.12: Closed tribotryoid

3.3.2 Sympodial inflorescences

Simple cymes (open)

In racemes, the apex of the main axis produces lateral branches andcontinues to grow. In contrast, the apex of the main axis in cymes turnsinto a flower shortly after a few lateral branches have been initiated.Their apices turn into flowers as well, and second-order branches takeover. In time, branches of higher and higher order are produced. Thus,the basic structure of a cymose inflorescence is captured by the partialL-system:

ω : ap1 : a → I[L]ap2 : a → I[L]Ap3 : A→ I[A]K

(3.5)

As in the open raceme, there is a single symbol with alternative ruleswhich specify that the vegetative apex a may change into a flower-producing apex A. Any one of the previously discussed mechanismsis available for timing this decision. Figure 3.13a shows an open cymewith branches curving in a spiral fashion, while Figure 3.13b shows one

3.3. Models of inflorescences 83

a b c

Figure 3.13: Open cymes: (a) spiral form, (b) zig-zag form, (c) double

with a zig-zag branching form.

Double cymes (open)

Frequently, not one but two lateral apices are produced under eachterminal apex as in Figure 3.13c. In this case the partial L-system is:

ω : ap1 : a → I[L]ap2 : a → I[L]Ap3 : A→ I[A][A]K

(3.6)

The two continuing lateral apices may develop at approximatelyequal rates (with the same plastochron) or with different rates, givingrise to asymmetric inflorescences. For example, the following L-systemscheme describes the development of rose campion (Lychnis coronaria) Lychnisas analyzed by Robinson [126]:

ω : A7

p1 : A7→ I[A0][A4]IK0

p2 : Ai→ Ai+1, 0 ≤ i < 7p3 : Ki→ Ki+1, i ≥ 0

Production p1 shows that at their creation time, the lateral apiceshave different states A0 and A4. Consequently, the first apex requireseight derivation steps to produce a flower and new branches, while thesecond requires only four steps. Each flower undergoes a sequence ofchanges, progressing from the bud stage to an open flower to a fruit.This developmental sequence is illustrated in Figure 3.14. According toproduction p1, the lateral apices branch at an angle of 45◦ and lie in aplane perpendicular to that defined by the mother axis and its sibling.Production p3 describes the linear elongation of internodes, while p4

84 Chapter 3. Developmental models

#include L(0),L(1),... /* leaf shapes */#include K(0),K(1),... /* flower shapes */

ω : A(7)p1 : A(t) : t=7→ FI(20)[&(60)∼L(0)]/(90)[&(45)A(0)]/(90)

[&(60)∼L(0)]/(90)[&(45)A(4)]FI(10)∼K(0)p2 : A(t) : t<7→ A(t+1)p3 : I(t) : t>0→ FFI(t-1)p4 : L(t) : * → L(t+1)p5 : K(t) : * → K(t+1)

Figure 3.14: Development of Lychnis coronaria

3.3. Models of inflorescences 85

a b c

Figure 3.15: Thyrsus: (a) spiral form, (b) zig-zag form, (c) double

and p5 capture the development of leaves and flowers over time. It isinteresting to note that at different developmental stages there are someopen flowers that have a relatively uniform distribution over the entireplant structure. This is advantageous to the plant since it increases thetime span over which seeds will be produced.

Cymes (closed)

Sympodial inflorescences that produce a terminal flower at some pointduring their development are called closed cymes. They result from theaddition of production

p4 : A → K

to the L-systems specified in (3.5) and (3.6), which define open singleand open double cymes.

Thyrsus (closed)

A thyrsus is an inflorescence with branches of cymes borne on a monopo-dially branching axis. Thus, it represents a mixed sympodial andmonopodial organization. Depending on the orientation of the flow-ers, a distinction between a thyrsus with cymes in a spiral form and ina zig-zag form can be made (Figure 3.15, a and b). Both of these typesare described by the following partial L-system:

86 Chapter 3. Developmental models

ω : ap1 : a → I[L]ap2 : a → I[L]Ap3 : A→ I[L][B]Ap4 : A→ Kp5 : B→ I[B]Kp6 : B→ K

In addition, a thyrsus may have double cymes (Figure 3.15c). In theclosed structure there are three developmental transformations. Thefirst represents the change from vegetative to flowering development onthe main axis (production p2). The second is necessary for the closureof the main axis with a terminal flower (p4). Both switches are relatedto the monopodial development of the main axis. The third transfor-mation is responsible for the formation of the flowers that terminatethe development of the sympodial structures (p6).

3.3.3 Polypodial inflorescences

Panicle

The term polypodial is not used in the botanical literature but is coinedhere to draw attention to the type of branching that represents con-tinuing development of the main axis as well as of the lateral apicesof a branch. The corresponding inflorescence type is usually called apanicle. The presence of two continuing apices at each new node isexpressed by the following production:

A → I[L][A]A

Since there can be nodes near the base of the plant that do not bearbranches, the usual initial rules are included to model the transitionfrom a purely vegetative to a flowering state. The resulting partialL-system is:

ω : ap1 : a → I[L]ap2 : a → I[L]Ap3 : A→ I[L][A]Ap4 : A→ K

An example of a paniculate structure is shown in Figure 3.16. Notethe presence of higher order branching and the lack of terminal racemes.Due to the repetitive application of production p3 at various levels ofbranching, the resulting structure is highly self-similar. The model in-cludes only two types of developmental transformations: the switch

3.3. Models of inflorescences 87

Figure 3.16: Panicle (elongated form)

from purely vegetative growth to the formation of the branching struc-ture (production p2), and the creation of terminal flowers (p4). Thetiming of the last production determines the flowering sequence of theplant. Two possible control mechanisms will be examined in detail, us-ing developmental models of the branching part of wall lettuce (Mycelismuralis) as examples.

The development of Mycelis is difficult to model for two reasons. MycelisFirst, the plant exhibits a basipetal flowering sequence, which meansthat flowering starts at the top of the plant and proceeds downwards.Secondly, at some developmental stages the plant has an acrotonicstructure, where the upper branches are more developed than the lowerones. Both phenomena are in a sense counter-intuitive, since it wouldseem that the older branches situated near the plant base should startgrowing and producing flowers before the younger ones at the planttop. To explain these effects, several models were proposed and for-mally analyzed by Janssen and Lindenmayer [77]. Their model II isrestated here as parametric L-system 3.2.

The axiom consists of three components. Modules F and A(0) rep- Model IIresent the initial segment and the apex of the main axis. Module I(20)is the source of a signal representing florigen. In nature, florigen issent towards the apex by leaves located at the plant base, which is notincluded in this model.

The developmental process consists of two phases that take placealong the main axis and are repeated recursively in branches of higherorders. First, the main axis is formed in a process of subapical growth

88 Chapter 3. Developmental models

#include O /* flower shape specification */#ignore / + ∼ O

ω : I(20)FA(0)p1 : S < A(t) : * → T(0)∼Op2 : A(t) : t>0→ A(t-1)p3 : A(t) : t=0→ [+(30)G]F/(180)A(2)p4 : S < F : * → FSp5 : F > T(c) : * → T(c+1)FU(c-1)p6 : U(c) < G : * → I(c)FA(2)p7 : I(c) : c>0→ I(c-1)p8 : I(c) : c=0→ Sp9 : S : * → εp10 : T(c) : * → ε

L-system 3.2: Mycelis muralis – Model II

specified by production p3. The apex produces consecutive segmentsF at the rate of one segment every three derivation steps (the delay iscontrolled by production p2), and initiates branches G positioned at anangle of 30◦ with respect to the main axis. The symbol G is interpretedhere in the same way as F . At this stage, the branches do not developfurther, which simulates the effect of apical dominance or the inhibitionof branch development during the active production of new branchesby the apex.

After a delay of 20 derivation steps, counted using production p7,an acropetal flower-inducing signal S is sent by production p8. Produc-tion p4 transports S across the segments at the rate of one internodeper step. Since new internodes are produced by the apex at a threetimes slower rate, the signal eventually reaches the apex. At this point,the second developmental phase begins. Production p1 transforms apexA(t) into a bud O. Further branch production is stopped and a signalT (c) is sent towards the base in order to enable the development oflateral branches. Parameter c is incremented by production p5 eachtime signal T (c) traverses an internode. Subsequently, production p6

introduces the value of parameter c into the corresponding branches,using module U(c) as a carrier. The successor of production p6 hasthe same format as the axiom, thus module I(c) determines the delaybetween the initiation of branch development and time signal S, sentto terminate further internode creation. This delay c is smallest forthe top branches and increases towards the plant base. Consequently,parameter c can be interpreted as the growth potential of the branches,allowing lower branches to grow longer than the higher ones. On theother hand, the development of the upper branches starts sooner, thusin some stages they will be more developed than the lower ones, andthe flowering sequence will progress downwards, corresponding to ob-

3.3. Models of inflorescences 89

#include K /* flower shape specification */#consider M S T V

ω : I(20)FA(0)p1 : S < A(t) : * → TV Kp2 : V < A(t) : * → TV Kp3 : A(t) : t>0 → A(t-1)p4 : A(t) : t=0 → M[+(30)G]F/(180)A(2)p5 : S < M : * → Sp6 : S > T : * → Tp7 : T < G : * → FA(2)p8 : V < M : * → Sp9 : T > V : * → Wp10: W : * → Vp11: I(t) : t>0 → I(t-1)p12: I(t) : t=0 → S

L-system 3.3: Mycelis muralis – Model III

servations of the real plant [77].A diagrammatic developmental sequence of Mycelis muralis simu-

lated using L-system 3.2 is shown in Figure 3.17. Initially, the segmentsare shown as bright green. The passage of florigen S turns them purple,and the lifting of apical dominance changes their color to dark green.Figure 3.18 represents a three-dimensional rendering of the same model.The three-dimensional structure differs from the two-dimensional di-agram only in details. The angle value associated with the module“/” in production p3 has been changed to 137.5◦, resulting in a spiralarrangement of lateral branches around the mother axis. The leavessubtending branches have been included in the model, and flowers havebeen assumed to undergo a series of changes from bud to open flowerto fruit.

Another developmental model of Mycelis, referred to here as model Model IIIIII, is given by L-system 3.3. The initial phases of development are thesame as in model II. First, apex A creates the main axis and initiateslateral branches (productions p3 and p4). Symbol M in the successor ofproduction p4 marks consecutive branching points. After a delay of 20steps (ω) counted by production p11, flowering signal S is generated atthe inflorescence base (p12) and sent up along the main axis (p5). Uponreaching the apex, S induces its transformation into a terminal flowerK, and initiates two basipetal signals T and V (p1). The basipetalsignals also can be initiated by production p2, which is needed for theproper timing of signals in the topmost lateral branch. Signal T prop-agates basipetally at the rate of one internode per derivation step (p6)and lifts apical dominance, thus allowing the lateral branches to grow(p7). The presence of the second basipetal signal V is the distinctive

90 Chapter 3. Developmental models

Figure 3.17: Development of Mycelis muralis

3.3. Models of inflorescences 91

Figure 3.18: A three-dimensional rendering of the Mycelis model

feature of model III. Its role is to enable the formation of flowers onthe lateral branches by generating the flowering signal S at their bases(p8). Since signal V propagates down the main axis at the rate of oneinternode per two derivation steps (p9, p10), the interval between thelifting of apical dominance by signal T and induction of flowering signalS by signal V increases linearly towards the inflorescence base. Thisallows the lower branches to grow longer than the upper ones, resultingin a structure that is more developed near the base than near the apexin later developmental stages.

This entire control process repeats recursively for each axis: its apexis transformed into a flower by signal S, the growth of lateral axes issuccessively enabled by signal T , and the second basipetal signal V issent to induce the flowering signal S in the next-order axes. Conse-quently, a basipetal flowering sequence is observed along all axes of thepanicle.

Model II controls the flowering on lateral branches using growth Biologicalrelevancepotential c accumulated by signal T on its way down, while model

III employs the time interval between signals T and V for the samepurpose. Since both models produce identical developmental sequences,it is not possible to decide which one is more faithful to nature withoutgathering additional data related to plant physiology. Nevertheless, themodels clearly indicate that the flowering sequence of Mycelis cannotbe explained simply in terms of two commonly recognized mechanisms,

92 Chapter 3. Developmental models

n=10, δ = 60◦

#include K /* flower shape specification */

ω : A∼Kp1 : A : * → [-/∼K][+/∼K]I(0)/(90)Ap2 : I(t) : !(t=2) → FI(t+1)p3 : I(t) : t=2 → I(t+1)[-FFA][+FFA]

Figure 3.19: Lilac inflorescences

3.3. Models of inflorescences 93

Figure 3.20: Geometry of a lilac inflorescence: (a) the decussate branchingpattern, (b) the inflorescence skeleton without flowers

the flowering signal and the lifting of apical dominance. Another factor,whether it is an accumulated delay or a third signal, is needed. Themathematical models bring forward evidence and assist in formulatingplausible hypotheses related to the control mechanisms that may beemployed by nature. The final answer will require further study of thereal plant.

The models of Mycelis employ relatively complicated control pro- Lilaccesses to explain the developmental sequence of a plant. On the otherhand, if only a static image of a panicle in a particular developmentalstage is needed, much simpler L-systems can be used. The L-systemthat generates the lilac inflorescences shown in Figure 3.19 is an ex-ample. Production p1 describes the subapical development of an axis.Production p2 models linear elongation of internodes in time and in-troduces a delay before p3 creates the lateral axes. The rotation of theapex by 90◦ (p1) results in a decussate branching pattern with consec-utive pairs of (n + 1)-order axes lying in the planes that pass throughthe n-order axis and are perpendicular to each other (Figure 3.20). Ascene including lilac inflorescences is shown in Figure 3.21.

3.3.4 Modified racemes

There are four frequently encountered types of inflorescences that aremorphological modifications of racemes. Their mature forms are of aspecial kind and need to be specified separately.

94 Chapter 3. Developmental models

Figure 3.21: The garden of L

3.3. Models of inflorescences 95

a b

Figure 3.22: Umbels: (a) simple, (b) compound

Umbel

An umbel is characterized by more than two internodes attached to asingle node, resulting in a typical umbrella-like shape. In a simple um-bel there are flowers at the ends of the lateral internodes (Figure 3.22a),while in compound umbels the branching pattern is repeated recursivelya certain number of times (Figure 3.22b). The partial L-system for asimple umbel is

ω : Ap1 : A → I[IK]n

and for a compound umbel of recursion depth two is

ω : Ap1 : A → I[IB]kBp2 : B → I[IC]lCp3 : C → I[IK]m

This type of inflorescence is commonly found in the family Umbelliferae.For example, Figure 3.23 presents a model of a wild carrot. Note that Wild carrotthe size of leaves decreases towards the top of the plant, producing aphase effect similar to that observed in simple racemes. In contrast, themost developed inflorescence is placed at the top of the plant, indicatingdevelopmental control by a hormone similar to that observed in mints(Figure 3.11).

Spike

An elongated raceme with closely packed flowers is called a spike. Manygrasses and sedges have this kind of inflorescence (Figure 3.24a). SeeFigure 4.17 (page 117) for a realistic model.

96 Chapter 3. Developmental models

Figure 3.23: Wild carrot

Spadix

A fleshy elongated raceme is called a spadix, and is frequently found inthe family Araceae (Figure 3.24b).

Capitulum

A fleshy spherical or disk-shaped raceme is called a capitulum or “head.”The sunflower head is an inflorescence of this kind, the oldest flowers be-ing at the margin and the youngest at the center (Figure 3.24c). Mem-bers of the family Compositae commonly have this type of structure.One characteristic feature is the spatial arrangement of components,such as flowers or seeds, which form early discernible spiral patterns. Adetailed description of these patterns is presented in the next chapter.

3.3. Models of inflorescences 97

a b c

Figure 3.24: Modified racemes: (a) spike, (b) spadix, (c) capitulum

Chapter 4

Phyllotaxis

The regular arrangement of lateral organs (leaves on a stem, scales ona cone axis, florets in a composite flower head) is an important aspectof plant form, known as phyllotaxis. The extensive literature generatedby biologists’ and mathematicians’ interest in phyllotaxis is reviewedby Erickson [36] and Jean [78]. The proposed models range widely frompurely geometric descriptions (for example, Coxeter [17]) to complexphysiological hypotheses tested by computer simulations (Hellendoornand Lindenmayer [59], Veen and Lindenmayer [151], Young [163]). Thischapter presents two models suitable for the synthesis of realistic imagesof flowers and fruits that exhibit spiral phyllotactic patterns.

Both models relate phyllotaxis to packing problems. The first oper-ates in a plane and was originally proposed by Vogel [154] to describethe structure of a sunflower head. A further detailed analysis was givenby Ridley [124, 125]. The second model reduces phyllotaxis to theproblem of packing circles on the surface of a cylinder. Its analysis waspresented by van Iterson [75] and reviewed extensively by Erickson [36].

The area of phyllotaxis is dominated by intriguing mathematicalrelationships. One of them is the “remarkable fact that the numbers ofspirals which can be traced through a phyllotactic pattern are predom-inantly integers of the Fibonacci sequence” [36, p. 54]. For example,Coxeter [17] notes that the pineapple displays eight rows of scales slop-ing to the left and thirteen rows sloping to the right. Furthermore, it isknown that the ratios of consecutive Fibonacci numbers Fk+1/Fk con-verge towards the golden mean τ = (

√5 + 1)/2. The Fibonacci angle

360◦τ−2, approximately equal to 137.5◦, is the key to the first modeldiscussed in this chapter.

100 Chapter 4. Phyllotaxis

nr n ∼ √

n 1 +

137 5 . °n 2 +

137 5 . °

#define a 137.5 /* divergence angle */#include D /* disk shape specification */

ω : A(0)p1 : A(n) : * → +(a)[f(n∧0.5)∼D]A(n+1)

Figure 4.1: Pattern of florets in a sunflower head, according to Vogel’s for-mula

4.1 The planar model

In order to describe the pattern of florets (or seeds) in a sunflower head,Vogel’s formulaVogel [154] proposed the formula

φ = n ∗ 137.5◦, r = c√

n, (4.1)

where:

• n is the ordering number of a floret, counting outward from thecenter. This is the reverse of floret age in a real plant.

• φ is the angle between a reference direction and the position vec-tor of the nth floret in a polar coordinate system originating atthe center of the capitulum. It follows that the divergence an-gle between the position vectors of any two successive florets isconstant, α = 137.5◦.

• r is the distance between the center of the capitulum and thecenter of the nth floret, given a constant scaling parameter c.

4.1. The planar model 101

a b c

Figure 4.2: Generating phyllotactic patterns on a disk. These three patternsdiffer only by the value of the divergence angle α, equal to (a) 137.3◦, (b)137.5◦ (the correct value), and (c) 137.6◦.

The distribution of florets described by formula (4.1) is shown in Fig-ure 4.1. The square-root relationship between the distance r and the Model

justificationfloret ordering number n has a simple geometric explanation. Assum-ing that all florets have the same size and are densely packed, the totalnumber of florets that fit inside a disc of radius r is proportional tothe disk area. Thus, the ordering number n of the most outwardlypositioned floret in the capitulum is proportional to r2, or r ∼ √

n.The divergence angle of 137.5◦ is much more difficult to explain.

Vogel [154] derives it using two assumptions.

• Each new floret is issued at a fixed angle α with respect to thepreceding floret.

• The position vector of each new floret fits into the largest existinggap between the position vectors of the older florets.

Ridley [125] does not object to these basic assumptions, but indi-cates that they are insufficient to explain the origin of the Fibonacciangle, and points to several arbitrary steps present in Vogel’s deriva-tion. He describes the main weakness as follows:

While it is reasonable to assume that the plant could con-tain genetic information determining the divergence angleto some extent, it is completely impossible for this alone tofix the divergence angle to the incredible accuracy occurringin nature, since natural variation in biological phenomenais normally rather wide. For example, for the 55- and 89-parastichies to be conspicuous, as occurs in most sunflowerheads, d must lie between 21

55and 34

89, a relative accuracy of

one part in 1869.

The critical role of the divergence angle α is illustrated in Figure 4.2.

102 Chapter 4. Phyllotaxis

Figure 4.3: Close-up of a daisy capitulum

Although a comprehensive justification of Vogel’s formula may requirefurther research, the model correctly describes the arrangement of flo-rets visible in actual capitula. The most prominent feature is two setsof spirals or parastichies, one turning clockwise, the other counterclock-Parastichieswise, which are composed of nearest neighboring florets. The numberof spirals in each set is always a member of the Fibonacci sequence; 21and 34 for a small capitulum, up to 89 and 144 or even 144 and 233 forlarge ones. For example, the capitulum of a daisy (Figure 4.3) exhibits34 clockwise spirals and 21 counterclockwise spirals, with directionsdetermined by following the spirals outwards from the capitulum cen-ter. In the image of a domestic sunflower capitulum (Figure 4.4), onecan discern 34 spirals running clockwise and 55 spirals running counter-clockwise. The number of perceived spirals depends on the capitulumsize expressed in terms of the number of component florets. If the fieldof attention is limited to a circle approximately 2/3 the size of the en-tire sunflower capitulum in Figure 4.4, the number of discernible spiralsbecomes 34 and 21.

Figure 4.4: Domestic sunflower head �

4.1. The planar model 103

#define S /* seed shape */#define R /* ray floret shape */#include M N O P /* petal shapes */

ω : A(0)p1 : A(n) : * → +(137.5)[f(n∧0.5)C(n)]A(n+1)p2 : C(n) : n <= 440 → ∼Sp3 : C(n) : 440 < n & n <= 565 → ∼Rp4 : C(n) : 565 < n & n <= 580 → ∼Mp5 : C(n) : 580 < n & n <= 595 → ∼Np6 : C(n) : 595 < n & n <= 610 → ∼Op7 : C(n) : 610 < n → ∼P

104 Chapter 4. Phyllotaxis

The dependence of the number of parastichies on the size of thefield of attention is yet another intriguing aspect of spiral phyllotaxis,as pointed out in the following excerpt from a letter by Alan Turing 1

quoted in [18]:

According to the theory I am working on now there is acontinuous advance from one pair of parastichy numbers toanother, during the growth of a single plant. . . . You willbe inclined to ask how one can move continuously from oneinteger to another. The reason is this — on any specimenthere are different ways in which the parastichy numbers canbe reckoned; some are more natural than others. During thegrowth of a plant the various parastichy numbers come intoprominence at different stages. One can also observe thephenomenon in space (instead of in time) on a sunflower.It is natural to count the outermost florets as say 21 + 34,but the inner ones might be counted as 8 + 13. . . . I don’tknow any really satisfactory account, though I hope to getone myself in about a year’s time.

A complete model of a flower head, suitable for realistic image syn-thesis, should contain several organs of various shapes. This is easilyachieved by associating different surfaces with specific ranges of the in-dex n. For example, consider the L-system that generates the sunflowerhead (Figure 4.4). The layout of components is specified by produc-Sunflower headtion p1, similar to that of the L-system in Figure 4.1. Productions p2

to p7 determine colors and shapes of components as a function of thederivation step number. The entire structure shown in Figure 4.4 wasgenerated in 630 steps. Alternatively, random selection of similar sur-faces could have been employed to prevent the excessive regularity ofthe resulting image.

Other extensions to the basic model consist of varying organ orien-tation in space and changing their altitude from the plane of the headas a function of n. For example, the sunflower plants included in Fig-ure 4.5 have flowers in four developmental stages: buds, young flowersstarting to open, open flowers and older flowers where the petals beginto droop. All flowers are generated using approximately the same num-ber of florets. The central florets are represented by the same surface ateach stage. The shape and orientation of surfaces representing petalsvary from one stage to another. The plants have been modeled as di-botryoids, with a single signal inducing a basipetal flowering sequence,as described in the previous chapter.

1To computer scientists, Alan Turing is best known as the inventor of the Turingmachine [146], which plays an essential role in defining the notion of an algorithm.However, biologists associate Turing’s name primarily with his 1952 paper, “Thechemical basis of morphogenesis” [147], which pioneered the use of mathematicalmodels in the study of pattern formation and advocated the application of comput-ers to simulate biological phenomena.

4.1. The planar model 105

Figure 4.5: Sunflower field

106 Chapter 4. Phyllotaxis

Figure 4.6: Zinnias

The zinnias (Figures 4.6 and 4.7) illustrate the effect of changing aOther examplespetal’s altitude, size and orientation as a function of n. The height atwhich a petal is placed decreases by a small amount as n increases. Thesize of each successive petal is incremented linearly. The orientation isalso adjusted linearly by a small angle increment. Thus, petals withsmall values of index n are placed more vertically, while petals withlarger indices n are more horizontal. Although the family Compositaeoffers the most examples of phyllotactic patterns, the same model canbe applied to synthesize images of other flowers, such as water-lilies(Figures 4.8 and 4.9) and roses (Figure 4.10).

Figure 4.7: Close-up of zinnias �

4.1. The planar model 107

108 Chapter 4. Phyllotaxis

Figure 4.8: Water-lily

Figure 4.9: Lily pond

4.2. The cylindrical model 109

Figure 4.10: Roses

4.2 The cylindrical model

The spiral patterns evident in elongated organs such as pine cones, Basic modelfir cones and pineapples, can be described by models that positioncomponents, in this case scales, on the surface of a cylinder. Van Iterson[75] divides phyllotactic patterns on cylinders into simple and conjugateones. In the case of a simple arrangement, all components lie on a singlegenerative helix. In contrast, conjugate patterns consist of two or moreinterleaved helices. This paper discusses simple phyllotactic patternsonly. They are generally characterized by the formula

φ = n ∗ α, r = const, H = h ∗ n, (4.2)

where:

• n is the ordering number of a scale, counting from the bottom ofthe cylinder;

• φ, r and H are the cylindrical coordinates of the nth scale;

• α is the divergence angle between two consecutive scales (as inthe planar case, it is assumed to be constant); and

• h is the vertical distance between two consecutive scales (mea-sured along the main axis of the cylinder).

110 Chapter 4. Phyllotaxis

A parametric L-system that generates the pattern described by for-ImplementationusingL-systems

mula (4.2) is given in Figure 4.11. The operation of this L-systemsimulates the natural process of subapical development characterizedby sequential production of consecutive modules by the top part of thegrowing plant or organ. The apex A produces internodes f(h) alongthe main axis of the modeled structure. Associated with each intern-ode is a disk ∼D placed at a distance r from the axis. This offset isachieved by moving the disk away from the axis using the module f(r),positioned at a right angle with respect to the axis by &(90). The spi-ral disk distribution is due to the module /(a), which rotates the apexaround its own axis by the divergence angle in each derivation step.

In the planar model, the constant divergence angle α = 137.5◦ isAnalysis ofmodel geometry found across a large variety of flower heads. The number of perceived

parastichies is determined by the capitulum size, and it changes as thedistance from the capitulum center increases. In contrast, a phyllotacticpattern on the surface of a cylinder is uniform along the entire cylinderlength. The number of evident parastichies depends on the values ofparameters α and h. The key problem, both from the viewpoint ofunderstanding the geometry of the pattern and applying it to generatesynthetic images, is to express the divergence angle α and the verticaldisplacement h as a function of the numbers of evident parastichies en-circling the cylinder in the clockwise and counterclockwise directions. Asolution to this problem was proposed by van Iterson [75] and reviewedby Erickson [36]. Our presentation closely follows that of Erickson.

The phyllotactic pattern can be explained in terms of circles packedon the surface of the cylinder. An evident parastichy consists of asequence of tangent circles, the ordering numbers of which form anarithmetic sequence with difference m. The number m is referred to asthe parastichy order. Thus, the circles on the cylinder surface may bearranged in two congruent 2-parastichies, five congruent 5-parastichies,and so on. The angular displacement between two consecutive circlesin an m-parastichy is denoted by δm. By definition, δm belongs to therange (−π, π] radians. The relation between the angular displacementδm and the divergence angle α is expressed by the equation

δm = mα − ∆m2π, (4.3)

where ∆m is an integer called the encyclic number. It is the numberof turns around the cylinder, rounded to the nearest integer, whichthe generative helix describes between two consecutive points of them-parastichy.

Usually, one can perceive two series of parastichies running in oppo-site directions (Figure 4.11). The second parastichy satisfies an equa-tion analogous to (4.3):

δn = nα − ∆n2π (4.4)

Consider the m- and n-parastichies starting at circle 0. In their pathsacross the cylinder, they will intersect again at circle mn. Assume

4.2. The cylindrical model 111

#define a 137.5281 /* divergence angle */#define h 35.3 /* vertical displacement */#define r 500 /* component offset from main axis */#include D /* disk shape specification */

ω : Ap1 : A : * → [&(90)f(r)∼D]f(h)/(a)A

Figure 4.11: Parastichies on the surface of a cylinder and on the unrolledcylinder. The L-system generates the cylindrical pattern.

that m and n are relatively prime; otherwise the phyllotactic patternwould have to contain several circles lying at the same height H and,contrary to the initial assumption, would not be simple. The circlemn is the first point of intersection between the m-parastichy and then-parastichy above circle 0. Consequently, the path from circle 0 to mnalong the m-parastichy, and back to 0 along the n-parastichy, encirclesthe cylinder exactly once. The section of m-parastichy between circles 0and mn consists of n+1 circles (including the endpoints), so the angulardistance between the circles 0 and mn is equal to nδm. Similarly, thedistance between circles 0 and mn, measured along the n-parastichy,can be expressed as mδn. As a result,

nδm − mδn = ±2π. (4.5)

112 Chapter 4. Phyllotaxis

mn

0 0

mnh

βγ

mδn nδm

Figure 4.12: An opposite parastichy triangle (as in Erickson [36, Fig. 3.8]).The base is formed by the circumference of the cylinder. The sides areformed by the parastichies.

The signs in equation (4.5) correspond to the assumption that the spi-rals encircle the cylinder in opposite directions; thus one of the valuesδ is positive and the other one is negative. Substituting the right sidesof equations (4.3) and (4.4) for δm and δn yields

n∆m − m∆n = ±1. (4.6)

To further analyze the pertinent geometric relationships, the cylin-der is cut along the vertical line passing through the center of circle 0and “unrolled” (Figure 4.11). The two parastichies and the circumfer-ence of the cylinder passing through point 0 form a triangle as shownin Figure 4.12. The perpendicular to the base from point mn dividesthis triangle into two right triangles. If d denotes the diameter of thecircles, then

(nδm)2 + (mnh)2 = (nd)2

and(mδn)2 + (mnh)2 = (md)2.

The above system of equations can be solved with respect to h and d:

h =√

(δ2m − δ2

n)/(n2 − m2) (4.7)

d =√

(n2δ2m − m2δ2

n)/(n2 − m2) (4.8)

or, taking into consideration equation (4.5),

d =√

2π(nδm + mδn)/(n2 − m2). (4.9)

4.2. The cylindrical model 113

a b

c d

0 01

23

45

67

89

1011

1213

1415

0 01 23 4 56 78 9 1011 12 1314 15

8

0 01 23 4 56 78 9 1011 12 1314 15

0 01 23 4 56 78 9 1011 12 1314 15

Figure 4.13: Patterns of tangent circles drawn on the surface of a cylinderas a function of circle diameter

The problem is to determine values of δm and δn. They are notsimply functions of parameters m and n. Figure 4.13 shows that, for agiven m and n, the values of δm and δn can be chosen from a certainrange, yielding parastichies of different steepness. In order to determinethis range, observe that at its limits the phyllotactic pattern changes;one previously evident parastichy disappears and another is formed.Thus, at the range limit, three evident parastichies coexist. It followsfrom Figure 4.13 that at one end of the range the third parastichy hasorder |m−n|, and at the other end it has order (m+n). Three coexistingparastichies imply that each circle is tangent to six other circles. Inother words, all circles lie in the vertices of a regular hexagonal grid,as seen in Figure 4.13a and c. Consequently, the angle β + γ at vertexmn (Figure 4.12) is equal to 2π/3. Expressing the base of the trianglein terms of its two sides and their included angle results in

(2π)2 = (nd)2 + (md)2 − 2(nd)(md) cos(2π/3)

or, after simplification,

d = 2π/√

m2 + mn + n2. (4.10)

Equations (4.9) and (4.10) yield

nδm + mδn = 2π(n2 − m2)/(m2 + mn + n2). (4.11)

114 Chapter 4. Phyllotaxis

Solving the system of equations (4.5) and (4.11) with respect to δm

and δn produces

δm = π(m + 2n)/(m2 + mn + n2) (4.12)

andδn = π(2m + n)/(m2 + mn + n2). (4.13)

Given the values of δm and δn, the divergence angle α can be foundfrom either equation (4.3) or (4.4), assuming that the encyclic numbers∆m or ∆n are known. It follows from the definition that these numbersare the smallest positive integers satisfying equation (4.6). A system-atic method for solving this equation, based on the theory of continuousfractions, is presented by van Iterson [75]. Erickson [36] points out thatin practice the solution can often be found by guessing. Another pos-sibility is to look for the smallest pair of numbers (∆m, ∆n) satisfying(4.6) using a simple computer program.

In conclusion, a phyllotactic pattern characterized by a pair of num-Patternconstruction bers (m,n) can be constructed as follows:

1. Find ∆m and ∆n from equation (4.6).

2. Find the range of admissible values of the angular displacementsδm and δn. The limits can be obtained from equations (4.12) and(4.13) using the values of m and n for one limit, and the pair(min{m,n}, |m − n|) for the other.

3. For a chosen pair of admissible displacement values δm and δn,calculate the divergence angle α from equation (4.3) or (4.4) andthe vertical displacement h from equation (4.8).

4. Find the diameter d of the circles from equation (4.8).

The diameter d does not enter directly in any formula used forimage synthesis, but serves as an estimate of the size of surfaces tobe incorporated in the model. This algorithm was applied to produceTable 4.1 showing parameter values for which three parastichies coexist.Triple-contact

patterns Given a pattern with two parastichies, this table provides the limits ofthe divergence angle α. For example, a (5,8) pattern can be formed forvalues of α ranging from 135.918365◦ to 138.139542◦, which correspondto the patterns (3,5,8) and (5,8,13), respectively.

Further information relating the divergence angle α to the verti-cal displacement h for various phyllotactic patterns is shown in Fig-ure 4.14. The arcs represent parameters of patterns with two paras-tichies (m,n). The branching points represent parameters of patternswith three parastichies (m,n,m + n).

4.2. The cylindrical model 115

m,n,m+n α (degrees) h d

(1, 1, 2) 180.000000 1.81380 -(1, 2, 3) 128.571426 0.777343 2.374821(1, 3, 4) 96.923073 0.418569 1.742642(2, 3, 5) 142.105270 0.286389 1.441462(1, 4, 5) 77.142860 0.259114 1.371104(3, 4, 7) 102.162163 0.147065 1.032949(3, 5, 8) 135.918365 0.111049 0.897598(2, 5, 7) 152.307693 0.139523 1.006115(1, 5, 6) 63.870968 0.175529 1.128493(4, 5, 9) 79.672134 0.089203 0.804479(4, 7, 11) 98.709671 0.058510 0.651536(3, 7, 10) 107.088600 0.068878 0.706914(3, 8, 11) 131.752579 0.056097 0.637961(5, 8, 13) 138.139542 0.042181 0.553204(5, 7, 12) 150.275223 0.049921 0.601820(2, 7, 9) 158.507462 0.081215 0.767613

Table 4.1: Cylinder formula values for triple-contact patterns

h

α0.00

0.10

0.20

0.30

0.40

0.50

0.60

0.70

0.80

0.90

1.00

1.10

1.20

1.30

1.40

1.50

1.60

1.70

1.80

60.0

65.0

70.0

75.0

80.0

85.0

90.0

95.0

100.

105.

110.

115.

120.

125.

130.

135.

140.

145.

150.

155.

160.

165.

170.

175.

180.

(1,2)

(1,2,3)

(1,3)

(1,2,3)

(2,3)

(1,3,4)(1,4) (1,3,4)

(3,4) (2,3,5)

(3,5)

(2,3,5)

(2,5)(1,4,5)(1,5) (1,4,5)

(4,5) (3,4,7)(4,7)

(3,4,7)

(3,7)(3,5,8)(3,8) (3,5,8)

(5,8)

(2,5,7)

(5,7)

(2,5,7)

(2,7)

Figure 4.14: The vertical displacement h as a function of the divergenceangle α for various phyllotactic patterns

116 Chapter 4. Phyllotaxis

Figure 4.15: Pineapples

Models of fruits synthesized using the cylindrical model are shown inFigures 4.15 and 4.16. The pineapple (Figure 4.15) is an examplePineappleof a pattern where a given scale has six neighbors, which belong to5-, 8- and 13-parastichies. The corresponding divergence angle α isequal to 138.139542◦. The spruce cones (Figure 4.16) were generatedSpruce conesusing the values m = 5, n = 8 and α = 137.5◦ (the divergence angleα for a (5, 8)-parastichy pattern belongs to the interval 135.918365◦ to138.139542◦). From these values, h and d were calculated as a functionof the radius of the cylinder. The effect of closing the bottom andtop of the pineapple and spruce cones was achieved by decreasing thediameter of the cylinder and the size of the scales.

A variant of the model of phyllotaxis on a cylinder can be used tomodel organs that are conical rather than cylindrical in shape. Forexample, Figure 4.17 shows a model of the sedge Carex laevigata. L-Sedgesystem 4.1 generates the male spike. Production p1 specifies the basic

4.2. The cylindrical model 117

Figure 4.16: Spruce cones

Figure 4.17: Carex laevigata: the male spike, the entire shoot, the femalespike

118 Chapter 4. Phyllotaxis

#define IRATE 1.025 /* internode growth rate */#define SRATE 1.02 /* spikelet growth rate */#include M /* spikelet shape specification */

ω : Ap1 : A : * → [&(5)f(1)∼M(1)]F(0.2)/(137.5)Ap2 : M(s) : s<3 → M(s*SRATE)p3 : f(s) : s<3 → f(s*SRATE)p4 : &(s) : s<15→ &(s*SRATE)p5 : F(i) : i<1 → F(i*IRATE)

L-system 4.1: The male spike of Carex laevigata

layout of the spikelets, similar to that given by the L-system in Fig-ure 4.11. The top portion of the spike has a conical shape due to thegrowth of spikelets for some time after their creation by the apex. Ac-cording to production p2, a spikelet grows by factor SRATE in eachderivation step, until it reaches the threshold size of 3. In an analogousway, productions p3, p4 and p5 capture the distance increase betweenspikelets and the spike axis, the increase of the branching angle, andthe elongation of internodes.

The models of phyllotaxis deal with the arrangement of organs inspace. For the purpose of mathematical analysis their shape is reducedto a simple geometric figure, usually a circle. However, in realisticimages the exact shape of the organs must be captured as well. Severaltechniques suitable for this purpose are outlined in the next chapter.

Chapter 5

Models of plant organs

Many concepts presented in the previous chapters were illustrated usingrealistic images, but the modeling techniques for leaves and petals havenot been described yet. Three approaches are discussed below.

5.1 Predefined surfaces

The standard computer graphics method for defining arbitrary surfaces Bicubic patchesmakes use of bicubic patches [9, 10, 40]. A patch is defined by threepolynomials of third degree with respect to parameters s and t. Thefollowing equation defines the x coordinate of a point on the patch:

x(s, t) = a11s3t3 + a12s

3t2 + a13s3t + a14s

3

+ a21s2t3 + a22s

2t2 + a23s2t + a24s

2

+ a31st3 + a32st

2 + a33st + a34s+ a41t

3 + a42t2 + a43t + a44

Analogous equations define y(s, t) and z(s, t). All coefficients are de-termined by interactively designing the desired shape on the screenof a graphics workstation. Complex surfaces are composed of severalpatches.

The surfaces are incorporated into a plant model in a manner sim- Turtleinterpretationilar to subfigures (Section 1.4.2). The L-system alphabet is extended

to include symbols representing different surfaces. When the turtleencounters such a symbol preceded by a tilde (∼), the correspondingsurface is drawn.

The exact position and orientation of surface S representing an ap-pendage is determined using a contact point PS, the heading vector �HS

and the up vector �V S as a reference (Figure 5.1). The surface is trans-lated in such a way that its contact point matches the current position

120 Chapter 5. Models of plant organs

Figure 5.1: Specification of an appendage

of the turtle, and is rotated to align its heading and up vectors with thecorresponding vectors of the turtle. If a surface represents an internalpart of the structure such as an internode, a distinction between theentry and exit contact points is made.

The majority of organs presented in this book have been modeledExamplesthis way. The petals of sunflowers, zinnias, water-lilies and roses shownin Chapter 4 provide good examples. Figure 5.2 illustrates an additionalimprovement in the appearance of organs, made possible by the appli-cation of textures to the surfaces of leaves, flowers and vine branches.

5.2 Developmental surface models

Predefined surfaces do not “grow.” String symbols can be applied tocontrol such features as the overall color and size of a surface, butthe underlying shape remains the same. In order to simulate plantdevelopment fully, it is necessary to provide a mechanism for changingthe shape as well as the size of surfaces in time. One approach is totrace surface boundaries using the turtle and fill the resulting polygons.Contour

tracing A sample L-system is given below:

ω : Lp1 : L → {−FX + X − FX − | − FX + X + FX}p2 : X → FX

Production p1 defines leaf L as a closed planar polygon. The braces{ and } indicate that this polygon should be filled. Production p2

increases the lengths of its edges linearly. The model of a fern shownFernin Figure 5.3 incorporates leaves generated using this method, withthe angle increment equal to 20◦. Note the phase effect due to the“growth” of polygons in time. A similar approach was taken to generatethe leaves, flowers and fruits of Capsella bursa-pastoris (Figure 3.5 onpage 74).

5.2. Developmental surface models 121

Figure 5.2: Maraldi figure by Greene [54]

Figure 5.3: The fern

122 Chapter 5. Models of plant organs

Figure 5.4: Surface specification using a tree structure as a framework

In practice, the tracing of polygon boundaries produces acceptable ef-Frameworkapproach fects only in the case of small, flat surfaces. In other cases it is more

convenient to use a tree structure as a framework. Polygon vertices arespecified by a sequence of turtle positions marked by the dot symbol(.). An example is given in Figure 5.4a. The letter G has been usedinstead of F to indicate that the segments enclosed between the bracesshould not be interpreted as the edges of the constructed polygon. Thenumbers correspond to the order of vertex specification by the turtle.

Figure 5.5 shows the development of a cordate leaf modeled usingCordate leafthis approach. The axiom contains symbols A and B, which initiate theleft-hand and right-hand sides of the blade. Each of the productions p1

and p2 creates a sequence of axes starting at the leaf base and graduallydiverging from the midrib. Production p3 increases the lengths of theaxes. The axes close to the midrib are the longest since they werecreated first. Thus, the shape of this leaf is yet another manifestationof the phase effect. The leaf blade is defined as a union of trianglesrather than a single polygon. Such triangulation is advantageous ifthe blade bends, for example due to tropism (Chapter 2). Figure 5.4bprovides an additional illustration of the model by magnifying the leftside of the leaf after four derivation steps.

5.2. Developmental surface models 123

ω : [A][B]p1 : A → [+A{.].C.}p2 : B → [-B{.].C.}p3 : C → GC

Figure 5.5: Developmental sequence of a cordate leaf generated using anL-system

The described method makes it possible to define a variety of leaves. Simple leavesTheir shapes depend strongly on the growth rates of segments. Forexample, a family of simple leaves and the corresponding parametricL-system are shown in Figure 5.6.

According to production p1, in each derivation step apex A(t) ex-tends the main leaf axis by segment G(LA,RA) and creates a pair oflateral apices B(t). New lateral segments are added by production p2.Parameter t, assigned to apices B by production p1, plays the role of“growth potential” of the branches. It is decremented in each derivationstep by a constant PD, and stops production of new lateral segmentsupon reaching 0. Segment elongation is captured by production p3.

For the purpose of analysis, it is convenient to divide a leaf bladeinto two areas. In the basal area, the number of lateral segments isdetermined by the initial value of growth potential t and constant PD.Since the initial value of t assigned to apices B increases towards theleaf apex, the consecutive branches contain more and more segments.On the other hand, branches in the apical area exist for too short atime to reach their limit length. Thus, while traversing the leaf fromthe base towards the apex, the actual number of segments in a branchfirst increases, then decreases. As a result of these opposite tendencies,the leaf reaches its maximum width near the central part of the blade.

124 Chapter 5. Models of plant organs

n=20, δ=60◦

#define LA 5 /* initial length - main segment */#define RA 1 /* growth rate - main segment */#define LB 1 /* initial length - lateral segment */#define RB 1 /* growth rate - lateral segment */#define PD 1 /* growth potential decrement */

ω : {.A(0)}p1 : A(t) : * → G(LA,RA)[-B(t).][A(t+1)][+B(t).]p2 : B(t) : t>0 → G(LB,RB)B(t-PD)p3 : G(s,r) : * → G(s*r,r)

Figure 5.6: A family of simple leaves generated using a parametric L-system

Figure LA RA LB RB PD

a 5 1.0 1.0 1.00 0.00b 5 1.0 1.0 1.00 1.00c 5 1.0 0.6 1.06 0.25d 5 1.2 10.0 1.00 0.50e 5 1.2 4.0 1.10 0.25f 5 1.1 1.0 1.20 1.00

Table 5.1: Values of constants used to generate simple leaves

5.2. Developmental surface models 125

Figure 5.7: A rose in a vase

Table 5.1 lists the values of constants corresponding to particular shapes. Shape controlFor a given derivation length, the exact position of the branch with thelargest number of segments is determined by PD. If PD is equal to 0,all lateral branches have an unlimited growth potential, and the basalpart of the leaf does not exist (Figure 5.6a). If PD equals 1, the basaland apical parts contain equal numbers of lateral branches (Figures 5.6b and f). Finer details of the leaf shape are determined by the growthrates. If the main axis segments and the lateral segments have the samegrowth rates (RA = RB), the edges of the apical part of the leaf arestraight (Figures 5.6 a and b). If RA is less than RB, the segmentsalong the main axis elongate at a slower rate than the lateral segments,resulting in a concave shape of the apical part (Figures 5.6 c and f). Inthe opposite case, with RA greater than RB, the apical part is convex(Figures 5.6 d and e). The curvature of the basal edges can be analyzedin a similar way.

126 Chapter 5. Models of plant organs

n=25, δ=60◦

#define LA 5 /* initial length - main segment */#define RA 1.15 /* growth rate - main segment */#define LB 1.3 /* initial length - lateral segment */#define RB 1.25 /* growth rate - lateral segment */#define LC 3 /* initial length - marginal notch */#define RC 1.19 /* growth rate - marginal notch */

ω : [{A(0,0).}][{A(0,1).}]p1 : A(t,d) : d=0 → .G(LA,RA).[+B(t)G(LC,RC,t).}]

[+B(t){.]A(t+1,d)p2 : A(t,d) : d=1 → .G(LA,RA).[-B(t)G(LC,RC,t).}]

[-B(t){.]A(t+1,d)p3 : B(t) : t>0 → G(LB,RB)B(t-1)p4 : G(s,r) : * → G(s*r,r)p5 : G(s,r,t) : t>1 → G(s*r,r,t-1)

Figure 5.8: A rose leaf

Figure 5.7 shows a long-stemmed rose with the leaves modeled accord-Rose leafing to Figure 5.8. The L-system combines the concepts explored inFigures 5.5 and 5.6. The axiom contains modules A(0, 0) and A(0, 1),which initiate the left-hand and right-hand side of the leaf. The de-velopment of the left side will be examined in detail. According toproduction p1, in each derivation step apex A(t, 0) extends the midribby internodes G(LA,RA) and creates two colinear apices B(t) pointingto the left. Further extension of the lateral axes is specified by produc-tion p3. The leaf blade is constructed as a sequence of trapezoids, withtwo vertices lying on the midrib and the other two vertices placed atthe endpoints of a pair of lateral axes formed in consecutive derivationsteps. The module G(LC,RC, t) introduces an offset responsible forthe formation of notches at the leaf margin. Production p4 describesthe elongation of internodes responsible for overall leaf shape, whileproduction p5 controls the size of the notches. The development of theright side of the blade proceeds in a similar manner, with production p2

5.2. Developmental surface models 127

Figure 5.9: Surface specification using stacked polygons

recreating the midrib and creating lateral apices pointing to the right.The bending of the midrib to the right is a result of tropism.

In the examples discussed so far, the turtle specifies the vertices of Nestedpolygonsone polygon, then moves on to the next. Further flexibility in surface

definition can be achieved by interleaving vertex specifications for dif-ferent polygons. The turtle interpretation of the braces is redefined inthe following way. A string containing nested braces is evaluated us-ing two data structures, an array of vertices representing the currentpolygon and a polygon stack. At the beginning of string interpretation,both structures are empty. The symbols {, } and . are then interpretedas follows:

{ Start a new polygon by pushing the current polygon onthe polygon stack and creating an empty current polygon.

. Append the new vertex to the current polygon.

} Draw the current polygon using the specified vertices,then pop a polygon from the stack and make it the currentpolygon.

An example of string interpretation involving nested braces is given inFigure 5.9.

The above technique was applied to construct the flowers of the Lily-of-the-valleylily-of-the-valley shown in Figure 3.4 (page 72), and magnified in Fig-

ure 5.10. A flower is represented by a polygon mesh consisting of fivesequences of trapezoids spread between pairs of curved lines that em-anate radially from the flower base. A single sequence is generated bythe following L-system:

ω : [X(36)A]/(72)[X(36)B]p1 : A : ∗ → [&GA{.].p2 : B : ∗ → B&.G.}p3 : X(a) : ∗ → X(a + 4.5)

128 Chapter 5. Models of plant organs

Figure 5.10: Structure of a lily-of-the-valley flower

Productions p1 and p2 create two adjacent framework lines and markpolygon vertices consistently with Figure 5.9. Production p3 controlsthe angle at which the framework lines leave the stalk at the flowerbase.

5.3 Models of compound leaves

So far, the discussion of organ models has focused on the definitionof surfaces. However, in the case of highly self-similar structures, theindividual surfaces become inconspicuous, and the expression of thegeometric relationships between younger and older parts of the struc-ture becomes the key issue. For example, Figure 5.11 shows compoundSymmetric

branching leaves often found in the family Umbelliferae. According to produc-tion p2, the apex A(0) creates two segments F (1) and a pair of lateralapices A(D) in each derivation step. Production p1 delays the develop-ment of the daughter branches by D steps with respect to the motherbranch. This pattern is repeated recursively in branches of higher or-der. Production p3 gradually elongates the internodes, and in this wayestablishes proportions between parts of a leaf. The values of the con-

Figure D R Derivation length

a 0 2.00 10b 1 1.50 16c 2 1.36 21d 4 1.23 30e 7 1.17 45

Table 5.2: Values of constants used to generate compound leaves

5.3. Models of compound leaves 129

#define D 1 /* apical delay */#define R 1.5 /* internode elongation rate */

ω : A(0)p1 : A(d) : d > 0 → A(d-1)p2 : A(d) : d = 0 → F(1)[+A(D)][-A(D)]F(1)A(0)p3 : F(a) : * → F(a*R)

Figure 5.11: Compound leaves

130 Chapter 5. Models of plant organs

#define D 1 /* apical delay */#define R 1.36 /* internode elongation rate */

ω : A(0)p1 : A(d) : d>0 → A(d-1)p2 : A(d) : d=0 → F(1)[+A(D)]F(1)B(0)p3 : B(d) : d>0 → B(d-1)p4 : B(d) : d=0 → F(1)[-B(D)]F(1)A(0)p5 : F(a) : * → F(a*R)

Figure 5.12: Compound leaves with alternating branching patterns

Figure D R Derivation length

a 1 1.36 20b 4 1.18 34c 7 1.13 46

Table 5.3: Values of the constants used to generate compound leaves withalternating branches

5.3. Models of compound leaves 131

stants used to generate the structures shown in Figure 5.11 are listedin Table 5.2. The model is sensitive to growth rate values — a changeof 0.01 visibly alters proportions. Leaves of the wild carrot, shown inFigure 3.23 (page 96), correspond closely to case Figure 5.11b.

Another type of compound leaf, with alternating branches, is ex- Alternatingbranchingamined in Figure 5.12. The values of the constants used in the L-

system are specified in Table 5.3. Further examples can be found inthe next chapter, devoted to the animation of developmental processes,and Chapter 8, which ties developmental models with fractals.

Chapter 6

Animation of plantdevelopment

The sequences of images used in Chapters 3 and 5 to illustrate the de- Motivationvelopment of inflorescences and compound leaves suggest the possibil-ity of using computer animation to visualize plant development. Froma practical perspective, computer animation offers several advantagesover traditional time-lapse photography.

• Photography is sensitive to imperfections in the underlying ex-periment. A disease or even a temporary wilting of a plant mayspoil months of filming.

• In nature, developmental processes are often masked by otherphenomena. For example, the growth of leaves can be difficultto capture because of large changes in leaf positions during theday. Similarly, positions of tree branches may be affected bywind. Computer animation makes it possible to abstract fromthese distracting effects.

• Animation can be used when time-lapse photography is imprac-tical because of the long development time of plants (e.g. yearsin the case of trees).

• In time-lapse photography, the initial position of the camera andthe light conditions must be established a priori, without know-ing the final shape of the plant. In computer animation all de-velopmental stages of the modeled plant are known in advance,allowing for optimal selection of the model view.

• Animation can be applied to visualize developmental mechanismsthat cannot be observed directly in real plants, such as the prop-agation of hormones and nutrients.

134 Chapter 6. Animation of plant development

• Animation offers an unprecedented means for visualizing the de-velopment of extinct plants on the basis of paleobotanical data.

The original formalism of L-systems provides a model of develop-Discretecharacter ofL-systems

ment that is discrete both in time and space. Discretization in timeimplies that the model states are known only at specific time inter-vals. Discretization in space means that the range of model states islimited to a finite number. Parametric L-systems remove the latter ef-fect by assigning continuous attributes to model components. However,the model states are still known only in discrete time intervals. Thispresents a problem in animation, where a smooth progression of thedeveloping forms is desirable.

This last statement requires further clarification. The very natureof animation is to produce the impression of continuous motion by dis-playing a sequence of still frames, captured at fixed time intervals. Whyis a continuous model of plant development needed if it will be used togenerate a fixed sequence of images in the final account? Wouldn’t itbe enough to retain the standard definition of L-systems and assumetime slices fine enough to produce the desired progression of forms?This approach, while feasible and useful, has three major drawbacks.

• A model can be constructed assuming longer or shorter time in-tervals, but once the choice has been made, the time step is a partof the model and cannot be changed easily. From the viewpointof computer animation it is preferable to control the time step bya single parameter, decoupled from the underlying L-system.

• The continuity criteria responsible for the smooth progression ofshapes during animation can be specified more easily in the con-tinuous time domain.

• It would be conceptually elegant to separate model development,defined in continuous time, from model observation, taking placein discrete intervals.

A developmental process is viewed as consisting of periods of contin-uous module expansion delimited by instantaneous module divisions.Special conditions are imposed to preserve the shape and growth ratesof the organism during these qualitative changes. An analogy can bedrawn to the theory of morphogenesis advanced by Thom [142], whoviewed shape formation as a piecewise continuous process with singu-larities called catastrophes.

Formally, development taking place in continuous time is capturedusing the formalism of timed DOL-systems. The key difference betweenthese L-systems and the types of L-systems considered so far lies inthe definition of the derivation function. In “ordinary” L-systems, thederivation length is expressed as the number of derivation steps. Intimed DOL-systems, the derivation length is associated with the totalelapsed time since the beginning of the observation.

6.1. Timed DOL-systems 135

6.1 Timed DOL-systems

Let V be an alphabet and R the set of positive real numbers (including0). The pair (a, τ) ∈ V × R is referred to as the timed letter a, andthe number τ is called the age of a. A sequence of timed letters, x =(a1, τ1) . . . (an, τn) ∈ (V × R)∗, is called a timed word over alphabet V .

A timed DOL-system (tDOL-system) is a triplet G = 〈V, ω, P 〉, Definitionwhere

• V is the alphabet of the L-system,

• ω ∈ (V ×R)+ is a nonempty timed word over V , called the initialword, and

• P ⊂ (V × R) × (V × R)∗ is a finite set of productions.

Instead of writing ((a, β), (b1, α1) . . . (bn, αn)) ∈ P , the notation(a, β) → (b1, α1) . . . (bn, αn) is used. The parameter β is referred toas the terminal age of the letter a, and each parameter αi is the initialage assigned to the letter bi by production P . The following assump-tions are made:

C1. For each letter a ∈ V there exists exactly one value β ∈ R suchthat (a, β) is the predecessor of a production in P .

C2. If (a, β) is a production predecessor and (a, α) is a timed let-ter that occurs in the successor of some production in P , thenβ > α.

According to these conditions, each letter has a uniquely definedterminal age. Furthermore, an initial age assigned to a letter by aproduction must be smaller than its terminal age, i.e., its lifetime (β −α) must be positive.

Let (a, β) → (b1, α1) . . . (bn, αn) be a production in P . A function DerivationD : ((V ×R)+ ×R) → (V ×R)∗ is called a derivation function if it hasthe following properties:

A1. D(((a1, τ1) . . . (an, τn)), t) = D((a1, τ1), t) . . .D((an, τn), t)

A2. D((a, τ), t) = (a, τ + t), if τ + t ≤ β

A3. D((a, τ), t) = D((b1, α1) . . . (bn, αn), t − (β − τ)), if τ + t > β

A derivation in a timed DOL-system is defined in terms of two typesof time variables. Global time t is common to the entire word underconsideration, while local age values τi are specific to each letter ai.

136 Chapter 6. Animation of plant development

Axiom A1 identifies t as the variable that synchronizes the entire de-velopment, and specifies that the lineages of all letters can be consideredindependently from each other (thus, no interaction between letters isassumed). With the progress of time t, each letter “grows older” untilits terminal age is reached (axiom A2). At this moment subdivisionoccurs and new letters are produced with initial age values specified bythe corresponding production (axiom A3). Condition C1 guaranteesthat the subdivision time is defined unambiguously, hence the devel-opment proceeds in a deterministic fashion. Condition C2 guaranteesthat, for any value of time t, the recursive references specified by axiomA3 will eventually end.

The above concepts are examined by formulating a timed DOL-Anabaenasystem that simulates the development of a vegetative part of the An-abaena catenula filament. Given the discrete model expressed by equa-tion (1.1) on page 5, the corresponding tDOL-system is as follows:

ω : (ar, 0)p1 : (ar, 1) → (al, 0)(br, 0)p2 : (al, 1) → (bl, 0)(ar, 0)p3 : (br, 1) → (ar, 0)p4 : (bl, 1) → (al, 0)

(6.1)

In accordance with the discrete model, it is assumed that all cells havethe same lifetime, equal to one time unit. The derivation tree is shownin Figure 6.1. The nodes of the tree indicate production applicationsspecified by axiom A3, and the triangular “arcs” represent the continu-ous aging processes described by axiom A2. The vertical scale indicatesglobal time. For example, at time t = 2.75 the filament consists of threecells, bl, ar and ar, whose current age is equal to 0.75.

According to the definition of time intervals corresponding to axiomsA2 and A3, a production is applied after the age τ + t exceeds theterminal age. Consequently, at division time the “old” cells still exist.For example, at time t = 2.0 the filament consists of two cells, al andbr, both of age τ = 1.

The above L-system can be simplified by considering cells of type bas young forms of the cells of type a. This is suggested by Figure 6.1where cells b simply precede cells a in time. The simplified L-systemhas two productions:

p1 : (ar, 2) → (al, 1)(ar, 0)p2 : (al, 2) → (al, 0)(ar, 1)

(6.2)

The corresponding derivation tree starting from cell (ar, 1) is shown inFigure 6.2. Note the similarity to the tree from the previous example.

Whether a natural developmental process or its mathematical modelModelobservation is considered, the choice of observation times and the act of observation

should not affect the process itself. In other terms, each derived wordshould depend only on the total elapsed time t, and not on the partition

6.1. Timed DOL-systems 137

Figure 6.1: Derivation tree representing the continuous development of An-abaena catenula described by the L-system is equation (6.1). Sections of thetriangles represent cell ages.

Figure 6.2: Derivation tree representing the continuous development of An-abaena catenula corresponding to the rules specified in equation (6.2)

138 Chapter 6. Animation of plant development

of t into intervals. The following theorem shows that timed DOL-systems satisfy this postulate.Theorem. Let G = 〈V, ω, P 〉 be a tDOL-system, and x ∈ (V ×R)+ bea timed word over V . For any values of ta, tb ≥ 0, the following holds:

D(D(x, ta), tb) = D(x, ta + tb)

Proof. Let us first consider the special case where the word x consistsof a single timed letter, (a0, τ0), and all productions in the set P takesingle letters into single letters. According to condition C1, there existsa unique sequence of productions from P such that:

(ai, βi) → (ai+1, αi+1), i = 0, 1, 2, . . .

Let (ak, τk) be the result of the derivation of duration ta that startsfrom (a0, τ0). According to axioms A2 and A3, and assuming thatta > β0 − τ0, this derivation can be represented in the form

D ((a0, τ0), ta)

= D((a1, α1), ta − (β0 − τ0))

= D((a2, α2), ta − (β0 − τ0) − (β1 − α1))

= · · ·= D((ak, αk), ta − (β0 − τ0) − (β1 − α1) − . . . − (βk−1 − αk−1))

= (ak, τk),

where

τk = αk + [ta − (β0 − τ0) −k−1∑i=1

(βi − αi)].

Since the sequence of recursive calls can be terminated only by an ap-plication of axiom A2, the index k and the age τk satisfy the inequality

αk < τk ≤ βk.

Due to condition C2, such an index k always exists and is unique.Let us now consider a derivation of duration tb > βk − τk that

starts from (ak, τk). Following the same reasoning, the result can berepresented as (am, τm), where

τm = αm + [tb − (βk − τk) −m−1∑

i=k+1

(βi − αi)]

andαm < τm ≤ βm.

By substituting the value of τk into the formula for τm, we obtain aftersimple transformations

τm = αm + [(ta + tb) − (β0 − τ0) −m−1∑i=1

(βi − αi)].

6.2. Selection of growth functions 139

Thus, the timed letter (am, τm) also results from the derivation of du-ration ta + tb starting with (a0, τ0):

(am, τm) = D(D((a0, τ0), ta), tb) = D((a0, τ0)x, ta + tb).

So far, we have considered only the case

ta > β0 − τ0, tb > βk − τk.

Two other cases, namely,

0 ≤ ta ≤ β0 − τ0, tb > βk − τk

andta > β0 − τ0, 0 ≤ tb ≤ βk − τk

can be considered in a similar way. The remaining case,

0 ≤ ta ≤ β0 − τ0, 0 ≤ tb ≤ βk − τk,

is a straightforward consequence of condition C2. This completes theproof of the special case. In general, a derivation that starts from aword (a1, τ1) . . . (an, τn) can be considered as n separate derivations,each starting from a single letter. This observation applies not only tothe initial word specified at time t = 0, but also to any intermediateword generated during the derivation. Consequently, any path in thederivation tree can be considered as a sequence of mappings that takessingle letters into single letters. Application of the previous reasoningseparately to every path concludes the proof. �

6.2 Selection of growth functions

Timed L-systems capture qualitative changes in model topology corre-sponding to cell (or, in general, module) divisions, and return the ageof each module as a function of the global time t. In order to completemodel definition, it is also necessary to specify the shape of each moduleas a function of its age. Potentially, such growth functions can be es-timated experimentally by observing real organisms [72, 73]. However,if detailed data is not available, growth functions can be selected froman appropriate class by choosing parameters so that the animation issmooth. This approach can be viewed as more than an ad hoc tech-nique for constructing acceptable animated sequences. In fact, Thompresents it as a general methodology for studying morphogenesis [142,page 4]:

The essence of the method to be described here consistsin supposing a priori the existence of a differential modelunderlying the process to be studied and, without knowingexplicitly what the model is, deducing from the single as-sumption of its existence conclusions relating to the natureof the singularities of the process.

140 Chapter 6. Animation of plant development

A technique for computing parameters of growth functions in thecase of nonbranching filaments and simple branching structures is givenbelow.

6.2.1 Development of nonbranching filaments

In a simple case of geometric interpretation of timed L-systems, symbolsrepresent cells that elongate during their lifetime and divide upon reach-ing terminal age. Model geometry does not change suddenly, whichmeans that

• the length of each cell is a continuous function of time, and

• the length of a cell before subdivision is equal to the sum of thelengths of the daughter cells.

The above continuity requirements are formalized in the context of aContinuityrequirements tDOL-system G = 〈V, ω, P 〉 as follows:

R1. Let [αmin, β] describe the life span of a timed letter a ∈ V . Theage αmin is the minimum of the initial age values assigned to aby the axiom ω or by some production in P . The terminal ageβ is determined by the predecessor of the production acting onsymbol a. The growth function g(a, τ), which specifies the lengthof cell a as a function of age τ , must be a continuous function ofparameter τ in the domain [αmin, β].

R2. For each production (a, β) → (b1, α1) . . . (bn, αn) in P the follow-ing equality holds:

g(a, β) =n∑

i=1

g(bi, αi) (6.3)

In practice, requirement R1 is used to select the class of growth func-tions under consideration, and the equations resulting from requirementR2 are used to determine the parameters in function definitions.

For example, in the case of the timed L-system specified in equa-Linear growthtion (6.2), requirement R2 takes the form

g(ar, 2) = g(al, 1) + g(ar, 0)g(al, 2) = g(al, 0) + g(ar, 1). (6.4)

Let us assume that the growth functions are linear functions of time:

g(al, τ) = Alτ + Bl

g(ar, τ) = Arτ + Br(6.5)

6.2. Selection of growth functions 141

Figure 6.3: Diagrammatic representation of the development of Anabaenacatenula, with (a) linear and (b) exponential growth of cells

By substituting equations (6.5) into (6.4), we obtain

2Ar + Br = (1Al + Bl) + (0Ar + Br) or 2Ar = Al + Bl

2Al + Bl = (0Al + Bl) + (1Ar + Br) or 2Al = Ar + Br.

The desired continuity of development is provided by all solutions ofthis system. They can be expressed in terms of coefficient c, whichrelates the growth rate of cells al to that of cells ar:

Al = cAr

Bl = Ar(2 − c)Br = Ar(2c − 1)

Figure 6.3a illustrates the developmental process considered for c =1. The diagram is obtained by plotting the cells in the filament ashorizontal line segments on the screen. Colors indicate cell type andage. The observation time t ranges from 1 (at the top) to 7 (at thebottom), with increment ∆t = 0.009.

The slopes of the side edges of the diagram represent growth ratesof the entire structure. Notice that they remain constant in the peri-ods between cell divisions, then change. This effect is disconcerting inanimation, since the rate of organism growth suddenly increases witheach cell division. In order to prevent this, it is necessary to extend re-quirements R1 and R2 to a higher order of continuity N . Specifically,equation (6.3) takes the form

g(k)(a, β) =n∑

i=1

g(k)(bi, αi) for k = 0, 1, . . . , N, (6.6)

142 Chapter 6. Animation of plant development

where g(k)(a, τ) is the kth derivative of the growth function g(a, τ) withrespect to age τ .

In the case of Anabaena, an attempt to achieve first-order continu-ity assuming linear growth functions yields an uninteresting solution,g(al, τ) = g(ar, τ) ≡ 0. Thus, more complex growth functions must beExponential

growth considered, such as an exponential function that can be used to approx-imate the initial phase of sigmoidal growth. Assume that the growthfunction has the form

g(al, τ) = g(ar, τ) = AeBτ . (6.7)

The objective is to find values of parameters A and B that satisfyequation (6.6) for k = 0, 1. By substituting equation (6.7) into (6.6),we obtain

ABke2B = ABkeB + ABk, (6.8)

which implies that zero-order continuity yields continuity of infiniteorder in this case. Solution of equation (6.8) for any value of k yields

B = ln1 +

√5

2≈ 0.4812. (6.9)

Parameter A is a scaling factor and can be chosen arbitrarily. Thecorresponding diagrammatic representation of development is shown inFigure 6.3b. The side edges of the diagram, representing the growthrates of the whole structure, are smooth exponential curves.

6.2.2 Development of branching structures

The notions of tDOL-system and growth function extend in a straight-forward way to L-systems with brackets. For example, the followingtDOL-system describes the recursive structure of the compound leavesanalyzed in Section 5.3.

ω : (a, 0)p1 : (a, 1) → (s, 0)[(b, 0)][(b, 0)](a, 0)p2 : (b, β) → (a, 0)

According to production p1, apex a produces an internode s, two lateralapices b and a younger apex a. Production p2 transforms the lateralapices b into a after a delay β. The daughter branches recursively repeatthe development of the mother branch.

Let us assume that the leaf development is first-order continuous,yielding the following equations for k = 0, 1:

g(k)(a, 1) = g(k)(s, 0) + g(k)(a, 0) (6.10)

g(k)(b, 0) = 0 (6.11)

g(k)(b, β) = g(k)(a, 0) (6.12)

6.2. Selection of growth functions 143

Equation (6.10) provides continuity along an existing axis upon theformation of a new internode s. Equation (6.11) specifies that a newlyformed lateral apex b has zero length and zero growth rate. Equa-tion (6.12) guarantees smooth transformation of apices b into apices a.

Let us further assume that apices a and internodes s expand expo-nentially,

g(a, τ) = AaeBaτ (6.13)

g(s, τ) = AseBsτ , (6.14)

for τ ∈ [0, 1]. The expansion of the lateral apices cannot be describedby an exponential function, since it would not satisfy equations (6.11).Consequently, a polynomial growth function g(b, τ) is chosen. Equa-tions (6.11) and (6.12) fix the function’s endpoints and the tangents atthe endpoints. Thus, in general, g(b, τ) must be a polynomial of degreethree or more, such as

g(b, τ) = Abτ3 + Bbτ

2 + Cbτ + Db. (6.15)

The system of equations (6.10) through (6.12) is solved using the initialsize Aa and the growth rate coefficient Ba as independent variables.By substituting (6.13) and (6.14) into equations (6.10) for k = 0, 1, weobtain

As = Aa(eBa − 1)

Bs = Ba.

Equations (6.11) and (6.15) yield

Cb = Db = 0.

Finally, substitution of (6.13) and (6.15) into (6.12) results in

Ab =Aa

β3(βBa − 2)

Bb =Aa

β2(3 − βBa).

Figure 6.4 shows a sequence of images produced by this model usingvalues β = Aa = 1 and Ba = 0.48. The branching angles are equal to45◦. The observation time t ranges from 6.9 to 7.7, with an incrementof 0.2. Note the gradual formation of lateral segments.

In the examples considered above, modules are represented asstraight lines, with growth functions controlling line lengths. Otherparameters, such as the branching angle, the diameter of segments andthe size of predefined surfaces, can be controlled in an analogous way.Generally, any developmental model captured by an OL-system withturtle interpretation can be converted into a tDOL-system and ani-mated.

144 Chapter 6. Animation of plant development

Figure 6.4: Developmental sequence of a branching structure modeled usinga tDOL-system

Chapter 7

Modeling of cellular layers

String L-systems, the first formalism considered in this book, are suit-able for the modeling of nonbranching filaments such as Anabaenacatenula. The introduction of brackets extends the class of modeledstructures to axial trees. However, many structures found in botanyhave a more complex topology, which can only be described by graphswith cycles. The developmental surface models presented in Section 5.2make it possible to specify a limited class of these graphs. This chapterdescribes a more general approach and applies it to simulate the devel-opment of single-layered cellular structures such as those found in ferngametophytes, animal embryos and plant epidermis. All structures con-sidered are of microscopic dimensions and relatively undifferentiated,yet the presented methods may bring us closer to the modeling of morecomplex patterns, such as the venation of leaves.

The modeling method consists of two stages. First, the topologyof the cell division patterns are expressed using the formalism of mapL-systems, which allows for the formation of cycles in a structure. Atthis stage the neighborhood relations between cells are established, butthe cell shapes remain unspecified. Next, cell geometry is modeledusing a dynamic method that takes into account the osmotic pressureinside the cells and the tension of cell walls. The development can beanimated by considering periods of continuous cell expansion, delimitedby instantaneous cell divisions.

7.1 Map L-systems

From a mathematical perspective, cellular layers can be represented Maps asmodels of celllayers

using a class of planar graphs with cycles, called maps [148]. Nakamuraet al. [102] characterize them as follows:

• A map is a finite set of regions. Each region is surrounded bya boundary consisting of a finite, circular sequence of edges thatmeet at vertices.

146 Chapter 7. Modeling of cellular layers

• Each edge has one or two vertices associated with it. The one-vertex case occurs when an edge forms a loop. The edges cannotcross without forming a vertex and there are no vertices withoutan associated edge.

• Every edge is a part of the boundary of a region.

• The set of edges is connected. Specifically, there are no islandswithin regions.

A map corresponds to a microscopic view of a cellular layer. Regionsrepresent cells, and edges represent cell walls perpendicular to the planeof view. The internal components of a cell are not considered.

The process of cell division can be described by map rewriting. ThisMap rewritingnotion is an extension of string rewriting as discussed in Section 1.2.In general, map-rewriting systems are categorized as sequential or par-allel, and can be region-controlled or edge-controlled [87]. Since severalcells may divide concurrently, a parallel rewriting system is needed.The second categorization has to do with the form of rewriting rules,which may express cell subdivisions in terms of region labels or edgelabels. Both approaches are suitable for biological modeling purposes[22]. This chapter focuses on the edge-controlled formalism of BinaryPropagating Map OL-system with markers, or mBPMOL-systems. ItmBPMOL-

systems was proposed by Nakamura, Lindenmayer and Aizawa [102] as a refine-ment of the basic concept of map L-systems introduced by Lindenmayerand Rozenberg [91]. The name is derived as follows. A map OL-systemis a parallel rewriting system that operates on maps and does not al-low for interaction between regions. In other words, regions are modi-fied irrespective of what happens to neighboring regions (a context-freemechanism). The system is binary because a region can split into atmost two daughter regions. It is propagating in the sense that the edgescannot be erased, thus regions (cells) cannot fuse or die. The markersspecify the positions of inserted edges that split the regions.

The choice of mBPMOL-systems as a modeling tool has two jus-tifications. First, they are more powerful than other interactionlessmap-rewriting systems described in the literature [19, 22, 23]. In addi-tion, markers have a biological counterpart in preprophase bands of mi-crotubules, which coincide with the attachment sites for division wallsformed during mitosis [55]. It should be noted, however, that doublewall systems, introduced by J. and H. B. Luck [93], may be relativelyeasier to specify [23].

An mBPMOL-system G consists of a finite alphabet of edge labelsDefinitionΣ, a starting map ω with labels from Σ, and a finite set of edge pro-ductions P . In general, the edges are directed, which is indicated by aleft or right arrow placed above the edge symbol. In some cases, theedge direction has no effect on the system operation. Such an edgeis called neutral and no arrow is placed above the symbol denoting it.Each production is of the form A → α, where the directed or neutral

7.1. Map L-systems 147

Figure 7.1: Examples of edge productions

edge A ∈ Σ is called the predecessor, and the string α, composed ofsymbols from Σ and special symbols [, ], +, and −, is called the succes-sor. The symbols outside square brackets specify the edge subdivisionpattern. Arrows can be placed above edge symbols to indicate whetherthe successor edges have directions consistent with, or opposite to, thepredecessor edge. Pairs of matching brackets [ and ] delimit markers,which specify possible attachment sites for region-dividing walls. Themarkers are viewed as short branches that can be connected to form acomplete wall. The strings inside brackets consist of two symbols. Thefirst symbol is either + or −, indicating whether the marker is placedto the left or to the right of the predecessor edge. The second symbolis the marker label, with or without an arrow. The left arrow indicatesthat the marker is directed towards the predecessor edge, and the rightarrow indicates that the marker is directed away from that edge. If noarrow is present, the marker is neutral.

For example, in the production→A → →

D←C [−←

E ]→BF , the directed prede- Production

syntaxcessor A splits into four edges D, C, B and F , and produces a marker E(Figure 7.1a). Successor edges D and B have the same direction as A,edge C has the opposite direction, and F is neutral. Marker E is placedto the right of A and is directed towards A. Note that this same pro-duction could be written as

←A → F

←B[+

←E ]

→C

←D (Figure 7.1b). As an example

of a production with a neutral predecessor, consider A→→B[−←

B]x[+←B]

←B. In

this case the result of production application does not depend on theassumed direction of the predecessor edge (Figure 7.1c).

A derivation step in an mBPMOL-system consists of two phases. Derivation

• Each edge in the map is replaced by successor edges and markersusing the corresponding edge production in P .

• Each region is scanned for matching markers.

148 Chapter 7. Modeling of cellular layers

ω : ABABp1 : A → B[-A][+A]Bp2 : B → A

Figure 7.2: Example of a map L-system. In the first step, a distinctionis made between the edge-rewriting phase and the connection of matchingmarkers.

Two markers are considered matching if:

• they appear in the same region,

• they have the same label, and

• one marker is directed away from its incident edge while the otheris directed towards the edge, or both markers are neutral.

If a match is found, the markers are joined to create a new edge thatwill split the region. The search for matching markers ends with thefirst match found, even though other markers entering the same regionmay also form a match. From the user’s perspective, the system be-haves nondeterministically since it chooses the pair of markers to beconnected. The unused markers are discarded.

The operation of mBPMOL-systems is illustrated in the followingExamplesexamples. The L-system shown in Figure 7.2 has two productions.Production p1 creates markers responsible for region division, whileproduction p2 introduces a delay needed to subdivide the regions alter-nately by horizontal and vertical edges.

The L-system shown in Figure 7.3 is a modified version of the pre-vious one. The only difference is the addition of an edge x, whichseparates the markers in the successor of production p1. This edgecreates a Z-shaped offset between the inserted edges A. Z-offsets andsymmetric S-offsets (Figure 7.4) can be observed in many biologicalstructures [22, 92].

Figure 7.5 illustrates the operation of an mBPMOL-system withdirected edges. Productions p1 and p3 create markers. Production p4

7.1. Map L-systems 149

ω : ABABp1 : A → B[-A]x[+A]Bp2 : B → A

Figure 7.3: Example of a map L-system

Figure 7.4: Offsets between four regions that result from the division of tworegions sharing a common edge: (a) Z-offset, (b) S-offset

ω :→A→B→C→D

p1 :→A → →

D[-→A]

→B

p2 :→B → →

B

p3 :→C → →

B[-←A]

→B

p4 :→D → →

C

Figure 7.5: Example of a map L-system with directed edges

150 Chapter 7. Modeling of cellular layers

transforms edge D into C, so that in each derivation step there is a pairof edges A and C to which productions p1 and p3 apply. Production p2

indicates that edges B do not undergo further changes.1 The resultingstructure is that of a clockwise spiral.

7.2 Graphical interpretation of maps

Maps are topological objects without inherent geometric properties. Inorder to visualize them, some method for assigning geometric interpre-tation must be applied. In its description, the biologically-motivatedterms cell and wall will be used instead of their mathematical counter-parts, region and edge.

Siero, Rozenberg and Lindenmayer [131] proposed a method that,Wallsubdivision in the simplest case, is expressed by the following rules:

• walls are represented by straight lines,

• the starting map is represented by a regular polygon, bounded bythe walls specified in the axiom,

• when a production subdivides a wall, all successor walls are ofequal length, and

• the position of a wall resulting from the union of two matchingmarkers is based on the position of these markers.

This wall subdivision method was used to draw Figures 7.2, 7.3 and7.5. However, in a biological context it creates cells with shapes thatare seldom observed in nature.

De Does and Lindenmayer [24] proposed a center of gravity methodCenter ofgravity method that produces more realistic shapes. According to this method each

interior vertex of the map is placed at the center of gravity of its neigh-bors. Such positioning of vertices has a sound biological justification;it minimizes hypothetical forces acting along cell walls, thus bringingthe entire structure to a state of minimum energy. However, if all ver-tices were positioned this way, the entire structure would collapse. Tocounteract this effect, the vertices on the map perimeter are displacedoutward a fixed distance. Unfortunately, this lacks biological justifi-cation and introduces sudden shape changes after cell divisions haveoccurred, making it unsuitable for animation purposes.

Assuming a dynamic point of view, the shape of cells and thus theDynamicmethod shape of the entire organism results from the action of forces. The

unbalanced forces due to cell divisions cause the gradual modificationof cell shapes until an equilibrium is reached. At this point, new celldivisions occur, and the search for an equilibrium resumes.

1In further L-systems such identity productions are omitted.

7.2. Graphical interpretation of maps 151

The dynamic method is based on the following assumptions:

• the modeled organism forms a single cell layer,

• the layer is represented as a two-dimensional network of massescorresponding to cell vertices, connected by springs that corre-spond to cell walls,

• the springs are always straight and obey Hooke’s law,

• each cell exerts pressure on its bounding walls directly propor-tional to the wall length and inversely proportional to the cellarea,

• the pressure on a wall is divided evenly between the wall vertices,

• the motion of masses is damped, and

• no other forces (for example, due to friction or gravity) are con-sidered.

The position of each vertex, and thus the shape of the layer, is com-puted as follows. As long as an equilibrium is not reached, unbalancedforces put masses in motion. The total force �F T acting on a vertex Xis given by the formula

�FT =∑

w∈W

�Fw + �Fd,

where �Fw are forces contributed by the set W of walls w incident toX, and �Fd = −b�v is a damping force, expressed as the product of adamping factor b and vertex velocity �v.

A wall w ∈ W contributes three forces acting on X (Figure 7.6).The tension �Fs acts along the wall, and its magnitude is determined byHooke’s law,

�Fs = −k(l − l0),

where k is the spring constant, l is the current spring length, and l0 isits rest length. The remaining forces, �P L and �P R, are due to the pressureexerted by the cells on the left and right sides of the wall. Each forceacts in a direction perpendicular to the wall and is distributed equallybetween the incident vertices. The magnitude of these forces equalsp · l, where p is the internal cell pressure and l is the wall length.The pressure is assumed to be inversely proportional to the cell area:p ∼ A−1. This assumption is derived from the equation describingosmotic pressure, p = SRT , as a function of the concentration of thesolute S (n moles per volume V ), the ideal gas constant R, and theabsolute temperature T [129]. Assuming that the cell volume V is

152 Chapter 7. Modeling of cellular layers

Figure 7.6: Forces acting on a cell vertex X

proportional to the area A captured by the two-dimensional modelunder consideration (V = Ah), pressure can be expressed as

p =nRT

Ah.

Thus p ∼ A−1, provided that the term nRT/h is constant. A convenientformula for calculating the area A is

A = |M∑i=1

(xi − xi+1)(yi + yi+1)/2|,

where (xi, yi) are coordinates of the M vertices of the cell, xM+1 = x1,and yM+1 = y1.

The force �F T acts on a mass placed at a map vertex. Newton’ssecond law of motion applies,

md2�x

dt2= �FT ,

where �x is the vertex position. If the entire structure has N vertices, asystem of 2N differential equations is obtained,

mid�vi

dt= �FTi

(�x1, · · · , �xN , �vi)

d�xi

dt= �vi,

where i = 1, 2, . . . , N . The task is to find the sequence of positions�x1, . . . , �xN at given time intervals, assuming that the functions �F Ti

andthe initial values of all variables �x 0

1 , . . . , �x 0N and �v 0

1 , . . . , �v 0N are known.

These initial values are determined as follows.

7.3. Microsorium linguaeforme 153

• Coordinates of the vertices of the starting map are included inthe input data for the simulation.

• Positions of existing vertices are preserved through a derivationstep. New vertices partition the divided walls into segments ofequal length. The initial velocities of all vertices are set to zero.

The system of differential equations with the initial values givenabove represents an initial value problem. It can be solved numericallyusing the forward (explicit) Euler method [41]. To this end, the dif-ferential equations are rewritten using finite increments ∆�vi, ∆�xi and∆t,

∆�v ki =

1

mi

�FTi

(�x k

1 , · · · , �x kN , �v k

i

)∆t

∆�x ki = �v k

i ∆t,

where the superscripts k = 0, 1, 2, . . . indicate the progress of time,t = k∆t. The position and velocity of a point i after time increment∆t are expressed as follows:

�v k+1i = �v k

i + ∆�v ki

�x k+1i = �x k

i + ∆�v kx

The iterative computation of the velocities �v ki and positions �x k

i is car-ried out for consecutive values of index k until all increments ∆�vi and∆�xi fall below a threshold value. This indicates that the equilibriumstate has been approximated to the desired accuracy. The next deriva-tion step is then performed. A system of equations corresponding tothe resulting map is created, and the search for an equilibrium stateresumes. In this way, the developmental process is simulated as periodsof continuous cell expansion, delimited by instantaneous cell divisions.Continuity of cell shapes during divisions is preserved by the rule thatsets the initial positions of vertices.

For example, Figure 7.7 illustrates the expansion of a structure gen-erated by the L-system specified in Figure 7.3 in a derivation of length4. Figure 7.7a shows the structure immediately after the insertion ofdivision walls, Figures 7.7b and c present intermediate wall positions,and Figure 7.7d describes the final structure at equilibrium.

7.3 Microsorium linguaeforme

In this section, the described simulation method is applied to visualizethe development of the fern gametophyte Microsorium linguaeforme.The biological data is based on observations conducted by de Boer [22].Fern gametophytes represent the sexually reproducing life stage of fern

154 Chapter 7. Modeling of cellular layers

a b

c d

Figure 7.7: Expansion of a cellular structure: (a) after division, (b) and (c)intermediate positions, (d) equilibrium

plants. They show no differentiation into stem, leaf and root, forminga plant body called a thallus. The development of a thallus can bedescribed conveniently in terms of the activity of the apical cell givingrise to segments, and the development of these segments. The modelingprocess captures repetitive patterns of cell divisions, so that large cel-lular structures can be described using a small number of productions.

The apical cell is the originator of the gametophyte structure. ItApical activitydivides repetitively, giving rise each time to a new apical cell and aprimary (initial) segment cell. The segment cells subsequently developinto multicellular segments. The division wall of an apical cell is at-tached to the thallus border on one side and to a previously createddivision wall on the other side. Thus, the division walls are oriented al-ternately to the left and to the right, yielding two columns of segmentsseparated by a zig-zag dividing line (Figure 7.8). The recursive natureof the apical activity can be expressed by the following formula (calleda cell division system [22]):

AL → SL | AR AR → AL | SR

This notation means that the cell on the left side of the arrow producestwo daughter cells separated by a wall.

7.3. Microsorium linguaeforme 155

Figure 7.8: Apical production of segments. The labels AR and AL denoteapical cells producing right segment SR and left segment SL, respectively.Dashed lines indicate the newly created division wall. The superscripts rep-resent segment age. The internal structure of segments is not shown.

Figure 7.9: Developmental sequence of a Microsorium segment

In describing the structure of a segment, we distinguish between pericli- Divisionpattern ofsegments

nal and anticlinal walls. Intuitively, periclinal walls are approximatelyparallel to the apical front of the thallus, and anticlinal walls are per-pendicular to this front. A more formal definition follows.

• In a primary segment, the apical front wall and one or more wallsopposing it are periclinal walls. The remaining walls are anticlinalwalls.

• A division wall attached to two periclinal walls is an anticlinalwall, and vice-versa.

In Microsorium, a wall is never attached to a periclinal wall on oneside and an anticlinal wall on the other side, so the above definitioncomprises all possible cases.

Microscopic observations of growing Microsorium gametophytes re-veal that most segments follow the same developmental sequence, showndiagrammatically in Figure 7.9. The primary segment cell Q1 is firstdivided by a periclinal wall into two cells, Q2 and Q3. Subsequently, the

156 Chapter 7. Modeling of cellular layers

C

C

C

C

Figure 7.10: Developmental sequence of a Microsorium gametophyte

7.3. Microsorium linguaeforme 157

ω :→A←Dx

→b

l1 :→a → ←

A [+←b ]

→i

l2 :→b → →

e [− →B ] x [+

→h ]

→d

l3 :→d → →

f

l4 :→f →

→g [− ←

h ] x [+→h ]

→d

l5 :→h → x [− →

f ] x

l6 :→i → →

c

l7 :→c → →

i [+←f ]

→i

l8 :→e → x [+ x] x

l9 :→g → x [− x] x [+ x] x

r1 :→A → ←

a [− ←B ]

→I

r2 :→B → →

E [+→b ] x [− →

H ]→D

r3 :→D → →

F

r4 :→F → →

G [+←H ] x [− →

H ]→D

r5 :→H → x [+

→F ] x

r6 :→I → →

C

r7 :→C → →

I [− ←F ]

→I

r8 :→E → x [− x] x

r9 :→G → x [+ x] x [− x] x

Figure 7.10 (continued): Developmental sequence of Microsorium

158 Chapter 7. Modeling of cellular layers

basal cell Q3 is divided by another periclinal wall into two “terminal”cells T that do not undergo further divisions. At the same time, the cellQ2 lying on the thallus border is divided by an anticlinal wall into twocells of type Q1. Each of these cells divides in the same way as the pri-mary cell. Consequently, the recursive nature of segment developmentcan be captured by the following cell division system:

Q1 →Q2

Q3

Q2 → Q1 | Q1 Q3 →T

T

In the above rules, a horizontal bar denotes a periclinal wall betweencells, and a vertical bar denotes an anticlinal wall.

The development of the Microsorium thallus is a result of concur-Development ofthe thallus rent divisions of the apical and segment cells. A single division of the

apical cell corresponds to a single step in segment development. A de-velopmental sequence that combines the activity of the apex and thesegments is shown in Figure 7.10. This figure also reveals offsets be-tween neighboring walls. It can be assumed that periclinal divisionwalls form S-offsets in the segments on the right side of the apex andZ-offsets in the segments on the left side.

In order to capture the development of Microsorium using the for-Map L-systemmalism of map L-systems, it is necessary to identify all combinationsof cells that may lie on both sides of a wall. Careful examinationof these combinations yields the wall labeling scheme shown in Fig-ure 7.10. Two walls have the same label if and only if they divide inthe same way.2 The uppercase letters apply to right segment walls, andthe corresponding lowercase letters denote symmetric walls in the leftsegments. A comparison of pairs of subsequent structures yields theL-system in the figure.

Apical cell division results from the application of productions r1

and l2 (creation of a right segment) or l1 and r2 (creation of a leftsegment). Subsequent cell divisions in right and left segments proceedsymmetrically. The development of a right segment is examined indetail.

The insertion of wall segment B creates the first cell Q1 of segmentSR

(1) (cf. Figure 7.8). Concurrently, wall D on the opposite side of thesegment is transformed into F (step 1). This transformation introducesa one-step delay before the application of production r4 which, togetherwith r2, splits cell Q1 into Q2 and Q3 by the first periclinal wall H (step2). As the derivation progresses, production r4 inserts subsequent per-iclinal walls H between pairs of anticlinal walls F (step 4). Production

2It is conceivable to formulate an algorithm that would assign labels consistentwith the above rule automatically. However, the labeling scheme given in Figure 7.10was obtained “by hand.”

7.3. Microsorium linguaeforme 159

Figure 7.11: Developmental sequence of a basal Microsorium segment

r3 introduces a delay needed to create walls F , which are inserted be-tween periclinal walls H and apical walls I using productions r5 and r7.Production r6 plays a role analogous to r3 — it introduces a one-stepdelay into the cycle which creates markers F at the apical front of thesegment. Thus, periclinal walls H and anticlinal walls F are producedalternately, in subsequent derivation steps. The last two productions,r8 and r9, create terminal walls x that do not undergo further changes.The first such wall is inserted between walls labeled G and E duringderivation step 3. The subsequent walls x are inserted every secondstep between pairs of walls G; only production r9 is applied in thesecases.

The L-system in Figure 7.10 was formulated under the assumption Basal segmentsthat all segments develop in the same way. However, in a real organismthe two oldest segments, situated at the thallus base, form a modifiedpattern with less extensive cell divisions. The developmental sequenceof a right basal segment is shown in Figure 7.11. The correspondingcell division system is given below.

Q1 →Q2

Q3

Q2 → Q1 | T Q3 →T

T

The development of a Microsorium gametophyte including basalsegments is captured by Figure 7.12. Only productions describing thedevelopment of the right side of the thallus are given. Their prede-cessors are denoted by uppercase letters. The corresponding lowercaseproductions, which complete the L-system definition, can be obtainedby switching the “case” of letters and the orientation of markers. Walldirections remain unchanged. For example, the right-side production

rx :→P →

←A[−

→b ]C

160 Chapter 7. Modeling of cellular layers

r1 :→A → ←

a [− ←B ]

→I

r2 :→B → →

E [+→b ] x [− →

H ]→D

r3 :→D → [− →

M ]→F

r4 :→F → →

G [+←H ] x [− →

H ]→D

r5 :→H → x [+

→F ] x

r6 :→I → →

C

r7 :→C → →

I [− ←F ]

→I

r8 :→E → x [− x] x

r9 :→G → x [+ x] x [− x] x

r10 :→J → →

L

r11 :→K → →

N

r12 :→L → x [+

←M ] x

r13 :→M → x [+

→L ] x

r14 :→N → →

O

r15 :→O → x [+

←L ]

→N

Figure 7.12: The initial map and productions for the right side of a Microso-rium linguaeforme with basal segments

corresponds to the left-side production

lx :→p → ←

a [+→B]c.

Additional examples can be found by comparing the left and rightcolumns of the L-system in Figure 7.10.

Assuming the starting map specified by Figure 7.12, a simulateddevelopmental sequence interpreted using the dynamic method to de-termine cell shape is given in Figure 7.13. Different colors are usedto indicate the apical cell, the alternating “regular” segments, and thebasal segments. A comparison of the last developmental stage with aModel

validation photograph of a real gametophyte (Figure 7.14) shows good correspon-dence between the model and reality with respect to structure topology,the relative sizes and shapes of cells, and the overall shape of the thal-lus. This result is particularly interesting from a biological perspective,since it indicates that genetically controlled cell division patterns playan important role in determining the shape of a structure.

7.3. Microsorium linguaeforme 161

Figure 7.13: Simulated development of Microsorium linguaeforme

Figure 7.14: Microphotograph of Microsorium linguaeforme at magnification70x

162 Chapter 7. Modeling of cellular layers

Figure 7.15: Developmental sequence of a Dryopteris segment

7.4 Dryopteris thelypteris

Gametophytes of other fern species follow a similar developmental pat-tern, with the apex producing segments alternately to the left andright. However, the cell division patterns within segments vary betweenspecies yielding different thalli shapes. For example, Figure 7.15 showssegment development in Dryopteris thelypteris. The corresponding celldivision system is given below.

Q1 →Q2

Q3

Q2 →Q4

Q3

Q3 →T

TQ4 → Q1 | Q1

A developmental cycle of length 3, starting at cell Q1, producestwo new cells Q1 separated by an anticlinal wall and a sequence offour terminal cells T separated by periclinal walls. A developmentalsequence that combines the activity of the apex and the segments isshown in Figure 7.17. As in the case of the map L-system in Figure 7.12,only productions for the right side of the thallus are shown.

Apical cell division results from the application of productions r1Map L-systemand l2 (creation of a right segment) or l1 and r2 (creation of a leftsegment). The subsequent cell divisions proceed in a symmetric wayin right and left segments. The development of a right segment isdescribed below.

The insertion of wall segment B creates the first cell Q1 of segmentSR

(1) (cf. Figure 7.8). Concurrently, wall D on the opposite side ofthe segment is transformed by production r6 into FG (step 1). Thistransformation introduces a one-step delay before the application ofproduction r7 which, together with r2, splits cell Q1 into Q2 and Q3 bythe first periclinal wall x (step 2). Meanwhile, production r3 replaceswall C by wall J (step 1), after which r4 replaces J by E (step 2). Thisintroduces a two-step delay before cell Q3 is subdivided periclinally intotwo cells T by production r5 (step 3). In the same step, productions r6

and r9 subdivide cell Q2 into cells Q3 and Q4, separated by periclinalwall H. Wall O from step 1 is transformed into wall R by productions

7.4. Dryopteris thelypteris 163

Figure 7.16: Simulated development of Dryopteris thelypteris

r15 (step 2) and r16 (step 3). Walls R and H are replaced by r17 andr14, resulting in the first anticlinal division of cell Q4 into two cells Q1

by wall I (step 4). At the same time, productions r5 and r6 split cellQ3 periclinally into two cells T . In the following derivation steps, eachof the newly created cells Q1 undergoes a sequence of changes similarto that described above. Production r8 introduces a one-step delaybefore Q1 is subdivided into Q2 and Q3 using r9 and r10 (analogousto r2 and r7). Productions r11 and r12 play a role similar to r5 andr6, while production r13 introduces a delay. Walls labeled x do notundergo further changes and cells T do not subdivide. A simulateddevelopmental sequence generated by the L-system in Figure 7.17 usingthe dynamic method to determine cell shape is given in Figure 7.16.

A comparison of Microsorium and Dryopteris gametophytes (either Microsoriumvs. Dryopterisreal or modeled) indicates that different division patterns of segment

cells have a large impact on the overall thalli shapes. In Microsorium,the number of marginal cells, situated at the apical front of a segment, isdoubled every second derivation step. The segments are approximatelytriangular, with a wide apical front, which results in the circular thallusshape. The apical front of Dryopteris segments is comparatively lessdeveloped. The number of marginal cells is doubled only every thirdstep, and the segments grow faster in length. The resulting thallusshape is concave.

164 Chapter 7. Modeling of cellular layers

(0)

b

C

D

RI

A

(3)

x

xx

x

x

x

x

x

x

x

x

x

Q12Q

3Q

T

T

a

IL

B

F

d

e

cx

O

R

H

G

JE

E

C

D3Q

Q4

p

e

k

j

(2)

RI

x

x

x x

Q1 2Q

3Q

x

A

o

P

b D

C

E

g

f

j

E

J

K

(1)

Q1

IL

a

O

Bd

c

e J

F

G

Figure 7.17: Developmental sequence of a Dryopteris gametophyte

7.4. Dryopteris thelypteris 165

ω :→A←D←C→b

r1 :→A → ←

a [− ←B ]

→O

r2 :→B → →

E [+→b ]

→C [− x]

→D

r3 :→C → ←

J

r4 :→J → →

E

r5 :→E → x [− x] x

r6 :→D → →

F [− →H ]

→G

r7 :→F → x [+ x] x [− x]

←J

r8 :→G → →

K

r9 :→K → →

E [+←H ]

→C [− x]

→D

r10 :→I → →

L [+ x] x [− x]→M

r11 :→L → x [+ x] x [− x] x

r12 :→M → →

L [+←H ] x [− →

H ]→N

r13 :→N → →

I

r14 :→H → x [+

→I ] x

r15 :→O → →

P

r16 :→P →

→Q

r17 :→Q → →

O [− ←I ]

→O

Figure 7.17 (continued): Developmental sequence of Dryopteris

166 Chapter 7. Modeling of cellular layers

In order to quantify the relationship between segment cell divisionpattern and thallus shape, de Boer [22] proposed an empirical measurecalled periclinal ratio. It is based on the following considerations:

• anticlinal growth takes place mainly along the margin and is ex-ponential since all marginal cells divide; and

• periclinal growth of a segment is linear, as cells displaced awayfrom the margin eventually stop dividing.

Since the division pattern is recursive, the average ratio of the numbersof marginal cells in neighboring segments converges to a constant A.Similarly, the difference between the numbers of cells along the pericli-nal boundary of two neighboring segments converges to a constant P .By computer simulation [22], it was found that if the periclinal ratioP/A is smaller than 1.25, the apical front is convex as in the case ofMicrosorium linguaeforme. In contrast, if P/A is larger than 2.0, thesimulated structure develops a concave apical front, corresponding to aheart-shaped thallus as in Dryopteris thelypteris. These results relatelocal growth to global shape. Periclinal ratios from 1.25 to 2.0 were notstudied in detail.

7.5 Modeling spherical cell layers

Although the scope of this book is limited to plants, it is interestingto note that the formalism of L-systems can be applied also to simu-late some developmental processes in animals. For example, during thecleavage stage of development, an animal embryo consists of a singleAnimal

embryos layer of cells that covers the surface of a spherical cavity. This structureis known as the blastula [6]. The cells divide synchronously in a regularpattern up to and including the 64-cell stage (6th cleavage). This devel-opment can be captured using an mBPMOL-system operating on thesurface of a sphere rather than on a plane. To this end, cell walls arerepresented as great circle arcs connecting vertices that are constrainedto the sphere surface.

The extension of the dynamic interpretation method from the planeto the surface of a sphere requires few changes. Osmotic pressure andwall tension are calculated as before. Since the resulting force maydisplace a vertex away from the surface of the sphere, the actual vertexposition is found by projecting the displaced point back to the sphere.During the cleavage stage, cells of embryos do not expand, thus theoverall size of the sphere is constant.

For example, deBoer [22] proposed the map L-system in Figure 7.18Patella vulgatato model the development of a snail embryo, Patella vulgata, accordingto data presented by Biggelaar [150]. The starting map and develop-mental sequence are shown in Figure 7.18, while Figure 7.19 presentsan alternative rendering. Each cell is approximated by a sphere cen-tered at the point obtained by raising the center of gravity of the cell

7.5. Modeling spherical cell layers 167

p1 : A → b[−a]x[+a]bp2 : a → B[+A]x[−A]Bp3 : B → ap4 : b → A

p5 :→C → →

D [+a]→E

p6 :→D → →

C [−A]x

p7 :→E → →

C

p8 : F → ←E [−a]G[+a]

→E

p9 : G → Jp10 : H → Ip11 : I → B[−A]x[+A]Bp12 : J → b[+a]x[−a]b

p13 : Z → →C [−F]H[+F]

←C

Figure 7.18: Developmental sequence of Patella vulgata (equatorial view)

168 Chapter 7. Modeling of cellular layers

vertices to the surface of the underlying spherical cavity. The radius isthe maximum distance from the center of gravity to the cell vertices.A comparison of the Patella model at the 16-cell stage (bottom leftof Figure 7.19) with an electron microscope image (Figure 7.20) showsgood correspondence between the model and reality.

7.6 Modeling 3D cellular structures

The previous sections presented a method for modeling cellular lay-Cellworksers extending in a plane or on the surface of a sphere. However, realcellular structures are three-dimensional objects. In order to capturethe three-dimensional aspect of cellular layers and model more complexstructures, Lindenmayer [85] proposed an extension of map L-systemscalled cellwork L-systems. The notion of a cellwork is characterized asfollows.

• A cellwork is a finite set of cells. Each cell is surrounded by oneor more walls (faces).

• Each wall is surrounded by a boundary consisting of a finite,circular sequence of edges that meet at vertices.

• Walls cannot intersect without forming an edge, although therecan be walls without edges (in the case of cells shaped as spheresor tori).

• Every wall is part of the boundary of a cell, and the set of wallsis connected.

• Each edge has one or two vertices associated with it. The edgescannot cross without forming a vertex, and there are no verticeswithout an associated edge.

• Every edge is a part of the boundary of a wall, and the set ofedges is connected.

Note that the terms cell and wall have different meanings in the contextof cellworks than in the context of maps.

The development of three-dimensional structures is captured usingmBPCOL-system an extension of mBPMOL-systems called marker Binary Propagating

Cellwork OL-systems [42]. An mBPCOL-system G is defined by a finitealphabet of edge labels Σ, a finite alphabet of wall labels Γ, a startingcellwork ω, and a finite set of edge productions P . The initial cellwork ωis specified by a list of walls and their bounding edges. As in the case ofmBPMOL-systems, edges may be directed or neutral. Each productionis of the form A : β → α, where the edge A ∈ Σ is the predecessor, thestring β ∈ {Γ+, ∗} is a list of applicable walls (* denotes all walls), andthe string α is the successor, which is composed of edge labels from Σ,

7.6. Modeling 3D cellular structures 169

Figure 7.19: Simulated development of Patella vulgata

Figure 7.20: An electron microscope image of Patella vulgata

170 Chapter 7. Modeling of cellular layers

Figure 7.21: The phases of a derivation step

wall labels from Γ and the symbols [ and ]. The symbols outside squarebrackets describe the subdivision pattern for the predecessor. Pairs ofmatching brackets [ and ] delimit markers that specify possible attach-ment sites for new edges and walls. As in the two-dimensional case,arrows indicate the relative directions of successor edges and markerswith respect to the predecessor edge. The list β contains all walls intowhich a marker should be inserted. In addition to the labels for edgesand markers, a production successor specifies the labels of walls thatmay be created as a result of a derivation step.

The syntax of a production is best explained using an example. TheProductionsyntax production

→A : 14 →

→D

←C2[

→E5]

→B3F

applies to an edge A if it belongs to one or more walls labeled 1 or4 (Figure 7.21a). The predecessor edge is subdivided into four edgesD, C, B and F . During a derivation step, marker E is introducedinto all walls of type 1 or 4 that share edge A (Figure 7.21b), andcan be connected with a matching marker inserted into the same wallby another production. As a result, the wall will split into two. Thedaughter wall having C as one of its edges will be labeled 2, and thewall having B as an edge will be labeled 3 (Figure 7.21c). Markers Ecan be connected only if both productions assign labels to the daughterwalls in a consistent way. Otherwise, the markers are considered non-matching and are discarded. If several walls bounding a cell split insuch a way that the sequence of new edges forms a closed contour, a

7.6. Modeling 3D cellular structures 171

Figure 7.22: Example of consistent edge productions

new wall bound by these edges may be created. In order for this tooccur, all markers involved must specify the same label for the newwall, 5 in this example (Figure 7.21d).

The limitation of the scope of a production to specific walls maycreate a consistency problem while rewriting edges. For instance, as-sume that walls 1 and 2 share edge A and the following productionsare in P :

p1 :→A : 1 →

→C

←E

p2 :→A : 2 → A

→B

Productions p1 and p2 are inconsistent since they specify two differentpartitions of the same edge. It is assumed that the mBPCOL-systemsunder consideration are free of such inconsistencies. This does notpreclude the possibility of applying several productions simultaneouslyto the same edge. For example, a production pair,

p1 :→A : 1 →

→C2[

→F 3]

←E4

p2 :→A : 2 →

→C5[

←D6]

←E7,

consistently divides edge A into segments C and E, although the mark-ers inserted into walls 1 and 2 are different (Figure 7.22).

According to the above discussion, a derivation step in an mBPCOL- Derivationsystem consists of three phases.

• Each edge in the cellwork is replaced by successor edges and mark-ers using one or more productions in P .

• Each wall is scanned for matching markers. If a match inducinga consistent labeling of daughter cells is found, the wall is subdi-vided. The selection of matching markers is done by the system.Unused markers are discarded.

• Each cell is scanned for a circular sequence of new division edges.If a cycle assigning the same label to the division wall is found, it isused to bound the wall that will divide the cell into two daughtercells. If different possibilities exist, the edges are selected by thesystem.

172 Chapter 7. Modeling of cellular layers

p1 : A : 1 → B1[A2]B1p2 : A : 2 → B2[C2]B2p3 : B : * → A

Figure 7.23: Example of a cellwork L-system

A wall may be subdivided more than once as long as new division edgesdo not intersect and a consistent labeling of daughter walls is possible.In contrast, a cell may be divided only once in any derivation step.

For example, Figure 7.23 presents a three-dimensional extension ofExamplethe map L-system from Figure 7.2. In the first derivation step, produc-tion p1 divides walls labeled 1, and production p2 divides walls labeled2. The inserted edges form a cycle that divides the cell with a new walllabeled 2. In the subsequent steps this process is repeated, generatinga pattern of alternating division walls. Production p3 introduces thenecessary delay.

The dynamic method for interpreting map L-systems is extended toDynamicinterpretation cellwork L-systems using the following assumptions:

• the structure is represented as a three-dimensional network ofmasses corresponding to cell vertices, connected by springs whichcorrespond to cell edges,

• the springs are always straight and obey Hooke’s law,

• for the purpose of force calculations, walls can be approximatedby flat polygons,

• the cells exert pressure on their bounding walls; the pressure ona wall is directly proportional to the wall area and inversely pro-portional to the cell volume,

7.6. Modeling 3D cellular structures 173

• the pressure on a wall is divided evenly between the wall vertices,

• the motion of masses is damped, and

• other forces are not considered.

The total force �F T acting on a vertex X is given by the formula

�FT =∑e∈E

�Fe +∑

w∈W

�Fw + �Fd,

where �Fe are forces contributed by the set of edges E incident to X,�Fw are forces contributed by the set of walls W incident to X, and�Fd = −b�v is a damping force. The forces �Fe act along the cell edges andrepresent wall tension. The forces �Fw are due to the pressure exertedby the cells on their bounding walls. The total force of pressure �Pexerted by a cell on a wall w has direction normal to w and is equalto p · A, where p is the internal cell pressure and A is the wall area.Calculation of the polygon area proceeds as in the two-dimensionalcase. The pressure p is assumed to be inversely proportional to thecell volume, p ∼ V −1, which corresponds to the equation describingosmotic pressure (Section 7.2). The volume V of a cell is calculated bytesselating the cell into tetrahedra. The resulting differential equationsare formed and solved as in the two-dimensional case.

A division pattern that frequently occurs in epidermal cell structures Epidermal cellsis described by the L-system in Figure 7.24, based on a cyclic cellworkL-system (a slightly different formalism) proposed by Lindenmayer [85].Productions p1, p2, p6 and p7 are responsible for cell divisions, while theremaining productions introduce delays such that the division patternis staggered.

On the surface, the cellular structures analyzed in this chapter mayappear quite unrelated to the models discussed previously. However,a closer inspection reveals many analogies. For example, consecutivelycreated segments of a fern gametophyte exhibit a phase effect corre-sponding to that observed in inflorescences. Furthermore, parts of anolder gametophyte situated near the apex have the same topology asthe entire thallus at an earlier developmental stage, which associatesthe recursive structures generated by map or cellwork L-systems withself-similar patterns created using string L-systems. As observed byOppenheimer [105], self-similarity appears to be one of the general prin-ciples organizing the world of botany. The next chapter discusses thistopic in more detail.

174 Chapter 7. Modeling of cellular layers

p1 : A : 123 → C3[E1]B2[D1]C3p2 : A : 4 → CB4[F1]C4p3 : B : ∗ → Ap4 : C : ∗ → Bp5 : E : ∗ → Dp6 : F : 123 → HGHp7 : F : 4 → H4[F1]G4[F1]H4p8 : G : ∗ → Fp9 : H : ∗ → G

Figure 7.24: Developmental sequence of epidermal cells: (a) The startingcellwork; (b), (d) and (f) cellworks immediately after cell divisions; (c), (e)and (g) the corresponding cellworks at equilibrium

Chapter 8

Fractal properties of plants

What is a fractal? In his 1982 book, Mandelbrot defines it as a set with Fractals vs.finite curvesHausdorff-Besicovitch dimension DH strictly exceeding the topological

dimension DT [95, page 15]. In this sense, none of the figures presentedin this book are fractals, since they all consist of a finite number ofprimitives (lines or polygons), and DH = DT . However, the situationchanges dramatically if the term “fractal” is used in a broader sense [95,page 39]:

Strictly speaking, the triangle, the Star of David, and thefinite Koch teragons are of dimension 1. However, bothintuitively and from the pragmatic point of view of the sim-plicity and naturalness of the corrective terms required, it isreasonable to consider an advanced Koch teragon as beingcloser to a curve of dimension log 4/log 3 than to a curveof dimension 1.

Thus, a finite curve can be considered an approximate renderingof an infinite fractal as long as the interesting properties of both areclosely related. In the case of plant models, this distinctive feature isself-similarity.

The use of approximate figures to illustrate abstract concepts has a Fractals vs.plantslong tradition in geometry. After all, even the primitives of Euclidean

geometry — a point and a line — cannot be drawn exactly. An in-teresting question, however, concerns the relationship between fractalsand real biological structures. The latter consist of a finite number ofcells, thus are not fractals in the strict sense of the word. To considerreal plants as approximations of “perfect” fractal structures would beacceptable only if we assumed Plato’s view of the supremacy of ideasover their mundane realization. A viable approach is the opposite one,to consider fractals as abstract descriptions of the real structures. Atfirst sight, this concept may seem strange. What can be gained by Complexity of

fractalsreducing an irregular contour of a compound leaf to an even more ir-regular fractal? Would it not be simpler to characterize the leaf using

176 Chapter 8. Fractal properties of plants

a smooth curve? The key to the answer lies in the meaning of the term“simple.” A smooth curve may seem intuitively simpler than a fractal,but as a matter of fact, the reverse is often true [95, page 41]. Accord-ing to Kolmogorov [80], the complexity of an object can be measuredby the length of the shortest algorithm that generates it. In this sense,many fractals are particularly simple objects.

The above discussion of the relationship between fractals and plantsPreviousviewpoints did not emerge in a vacuum. Mandelbrot [95] gives examples of the re-

cursive branching structures of trees and flowers, analyzes theirHausdorff-Besicovitch dimension and writes inconclusively “trees maybe called fractals in part.” Smith [136] recognizes similarities betweenalgorithms yielding Koch curves and branching plant-like structures,but does not qualify plant models as fractals. These structures are pro-duced in a finite number of steps and consist of a finite number of linesegments, while the “notion of fractal is defined only in the limit.” Op-penheimer [105] uses the term “fractal” more freely, exchanging it withself-similarity, and comments: “The geometric notion of self-similaritybecame a paradigm for structure in the natural world. Nowhere is thisprinciple more evident than in the world of botany.” The approach pre-sented in this chapter, which considers fractals as simplified abstractrepresentations of real plant structures, seems to reconcile these previ-ous opinions.

But why are we concerned with this problem at all? Does the no-Fractals inbotany tion of fractals provide any real assistance in the analysis and modeling

of real botanical structures? On the conceptual level, the distinctivefeature of the fractal approach to plant analysis is the emphasis onself-similarity. It offers a key to the understanding of complex-looking,compound structures, and suggests the recursive developmental mech-anisms through which these structures could have been created. Thereference to similarities in living structures plays a role analogous to thereference to symmetry in physics, where a strong link between conser-vation laws and the invariance under various symmetry operations canbe observed. Weyl [159, page 145] advocates the search for symmetryas a cognitive tool:

Whenever you have to deal with a structure-endowed entityΣ, try to determine its group of automorphisms, the groupof those element-wise transformations which leave all struc-tural relations undisturbed. You can expect to gain a deepinsight into the constitution of Σ in this way.

The relationship between symmetry and self-similarity is discussedin Section 8.1. Technically, the recognition of self-similar features ofplant structures makes it possible to render them using algorithms de-veloped for fractals as discussed in Section 8.2.

8.1. Symmetry and self-similarity 177

Figure 8.1: The Sierpinski gasket is closed with respect to transformationsT1, T2 and T3 (a), but it is not closed with respect to the set including theinverse transformations (b).

8.1 Symmetry and self-similarity

The notion of symmetry is generally defined as the invariance of a con-figuration of elements under a group of automorphic transformations.Commonly considered transformations are congruences, which can beobtained by composing rotations, reflections and translations. Couldwe extend this list of transformations to similarities, and consider self-similarity as a special case of symmetry involving scaling operations?

On the surface, this seems possible. For example, Weyl [159, page 68]suggests: “In dealing with potentially infinite patterns like band orna-ments or with infinite groups, the operation under which a pattern isinvariant is not of necessity a congruence but could be a similarity.”The spiral shapes of the shells Turritella duplicata and Nautilus aregiven as examples. However, all similarities involved have the samefixed point. The situation changes dramatically when similarities withdifferent fixed points are considered. For example, the Sierpinski gasketis mapped onto itself by a set of three contractions T1, T2 and T3 (Fig-ure 8.1a). Each contraction takes the entire figure into one of its threemain components. Thus, if A is an arbitrary point of the gasket, andT = Ti1Ti2 . . . Tin is an arbitrary composition of transformations T1, T2

and T3, the image T (A) will belong to the set A. On the other hand,if the inverses of transformations T1, T2 and T3 can also be includedin the composition, one obtains points that do not belong to the set Anor its infinite extension (Figure 8.1b). This indicates that the set oftransformations that maps A into itself forms a semigroup generatedby T1, T2 and T3, but does not form a group. Thus, self-similarity is aweaker property than symmetry, yet it still provides a valuable insightinto the relationships between the elements of a structure.

178 Chapter 8. Fractal properties of plants

Figure 8.2: The fern leaf from Barnsley’s model [7]

8.2 Plant models and iterated function sys-

tems

Barnsley [7, pages 101–104] presents a model of a fern leaf (Figure 8.2),generated using an iterated function system, or IFS. This raises a ques-tion regarding the relationship between developmental plant modelsexpressed using L-systems and plant-like structures captured by IFSes.This section briefly describes IFSes and introduces a method for con-structing those which approximate structures generated by a certaintype of parametric L-system. The restrictions of this method are ana-lyzed, shedding light on the role of IFSes in the modeling of biologicalstructures.

By definition [74], a planar iterated function system is a finite setIFS definitionof contractive affine mappings T = {T1, T2, . . . , Tn} which map theplane R×R into itself. The set defined by T is the smallest nonemptyset A, closed in the topological sense, such that the image y of anypoint x ∈ A under any of the mappings Ti ∈ T also belongs to A.It can be shown that such a set always exists and is unique [74] (seealso [118] for an elementary presentation of the proof). Thus, startingfrom an arbitrary point x ∈ A, one can approximate A as a set ofimages of x under compositions of the transformations from T . On

8.2. Plant models and iterated function systems 179

Figure 8.3: A comparison of three attracting methods for the rendering ofa set defined by an IFS: (a) deterministic method using a balanced tree ofdepth n = 9 with the total number of points N1 = 349, 525, (b) deterministicmethod using a non-balanced tree with N2 = 198, 541 points, (c) stochasticmethod with N3 = N2 points

the other hand, if the starting point x does not belong to A, the con-secutive images of x gradually approach A, since all mappings Ti arecontractions. For this reason, the set A is called the attractor of theIFS T . The methods for rendering it are based on finding the images Rendering

methodsTik(Tik−1(. . . (Ti1(x)) . . .)) = xTi1 . . . Tik−1

Tik , and are termed attractingmethods. According to the deterministic approach [123], a tree of trans-formations is constructed, with each node representing a point in A.Various strategies, such as breadth-first or depth-first, can be devised totraverse this tree and produce different intermediate results [60]. If thetransformations in T do not have the same scaling factors (Lipschitzconstants), the use of a balanced tree yields a non-uniform distribu-tion of points in A. This effect can be eliminated by constructing anon-balanced tree, using a proper criterion for stopping the extensionof a branch [60]. An alternative approach for approximating the set Ais termed the chaos game [7] (see also [107, Chapter 5]). In this case,only one sequence of transformations is constructed, corresponding toa single path in the potentially infinite tree of transformations. Thetransformation applied in each derivation step is selected at random.In order to achieve a uniform distribution of points in the attractor,the probability of choosing transformation Ti ∈ T is set according toits Lipschitz constant. Figure 8.3 illustrates the difference between thestochastic and deterministic methods of rendering the attractor. The

180 Chapter 8. Fractal properties of plants

underlying IFS consists of four transformations, given below using ho-mogeneous coordinates [40]:

T1 =

0.00 0.00 0.00

0.00 0.16 0.000.00 0.00 1.00

T2 =

0.20 0.23 0.00−0.26 0.22 0.00

0.00 1.60 1.00

T3 =

−0.15 0.26 0.00

0.28 0.24 0.000.00 0.44 1.00

T4 =

0.85 −0.04 0.00

0.04 0.85 0.000.00 1.60 1.00

Other methods for the rendering of the set A, defined by an inter-ated function system T , include the repelling or escape-time methodand the distance method [60, 118]. Both methods assign values topoints outside of A. The first method determines how fast a point isrepelled from A to infinity by the set of inverse transformations T−1

i ,where Ti ∈ T . An example of the application of this method, withescape time values represented as a height field, is shown in Figure 8.4.The second method computes the Euclidean distance of a point fromthe attractor A.

The problem of constructing an IFS that will approximate a branch-IFSconstruction ing structure modeled using an L-system can now be considered. This

discussion focuses specifically on structures that develop in a biologi-cally justifiable way, by subapical branching (Section 3.2). The com-pound leaf shown in Figure 5.11a on page 129 will be used as a workingexample. In this case, the apical delay D is equal to zero, and theL-system can be represented in the simplified form:

ω : Ap1 : A : ∗ → F (1)[+A][−A]F (1)Ap2 : F (a) : ∗ → F (a ∗ R)

(8.1)

Figure 8.4: Fern dune �

8.2. Plant models and iterated function systems 181

182 Chapter 8. Fractal properties of plants

Figure 8.5: Initial sequences of structures generated by the L-systems spec-ified in equations (8.1) and (8.2)

This L-system operates by creating segments of constant size, thenincreasing their length by constant factor R in each derivation step(Figure 8.5a). As discussed in Section 1.10.3, a structure with thesame proportions can be obtained by successively appending segmentsof decreasing length (Figure 8.5b):

ω : A(1)p1 : A(s) : ∗ → F (s)[+A(s/R)][−A(s/R)]F (s)A(s/R) (8.2)

Let An(s) denote the structure generated by module A(s) in n ≥ 1derivation steps. According to production p1, the following equalityholds:

An(s) = F (s)[+An−1(s/R)][−An−1(s/R)]F (s)An−1(s/R) (8.3)

8.2. Plant models and iterated function systems 183

Figure 8.6: Illustration of equation (8.4), with µ1 = F (1) − (45) and µ2 =[−(45)F (0.5)]F (0.5)

It is important to clearly distinguish between a parametric word µ (a Properties ofturtleinterpretation

string of modules) and its turtle interpretation J (µ) (a set of points inthe plane). The symbol M(µ) will be used to denote the transformationinduced by µ. This transformation moves the turtle from its initialposition and orientation to those resulting from the interpretation ofword µ. According to the definition of turtle interpretation (Chapter 1),if a word µ is decomposed into subwords µ1 and µ2 such that µ2 doesnot contain unbalanced right brackets, then

J (µ) = J (µ1µ2) = J (µ1) ∪ J (µ2)M(µ1). (8.4)

See Figure 8.6 for an illustration. By applying equation (8.4) to (8.3),we obtain

J (An(s)) = J (F (s)) ∪J (An−1(s/R))M(F (s)+) ∪J (An−1(s/R))M(F (s)−) ∪ (8.5)

J (F (s))M(F (s)) ∪J (An−1(s/R))M(F (s)F (s)),

which is true for any n ≥ 1. Now let A(s) be the limit of the sequence Passage toinfinityof sets J (An(s)) from equation 8.6,

A(s) = limn→∞J (An(s)).

At the limit we obtain:

A(s) = J (F (2s)) ∪ A(s/R)(T ′1 ∪ T ′

2 ∪ T ′3), (8.6)

184 Chapter 8. Fractal properties of plants

where T ′1 = T (F (s)+), T ′

2 = T (F (s)−), and T ′3 = T (F (2s)). Let

S(s/R) be the operation of scaling by s/R, then

A(s/R) = A(s)S(s/R).

By noting Ti = T ′iS(s/R) for i = 1, 2, 3, equation (8.6) can be trans-

formed to

A(s) = J (F (2s)) ∪ A(s)(T1 ∪ T2 ∪ T3). (8.7)

The solution of this equation with respect to A(s) is

A(s) = J (F (2s))(T1 ∪ T2 ∪ T3)∗, (8.8)

where (T1 ∪ T2 ∪ T3)∗ stands for the iteration of the union of transfor-

mations T1, T2 and T3. Equation (8.8) suggests the following methodfor constructing the set (A(s)):

• create segment J (F (2s))

• create images of J (F (2s)) using transformations T1, T2, T3 andtheir compositions

Equation (8.7) and the method of constructing the set A(s) basedon equation (8.8) are closely related to the definition of iterated functionsystems stated at the beginning of this chapter. However, instead ofstarting from an arbitrary point x ∈ A(s), the iteration begins withthe set J (F (2s)). Although this is simply a straight line segment, aquestion arises as to how its generation can be incorporated into anIFS. Two approaches can be distinguished.

The first approach is related to the notions of hierarchical iteratedControlled IFSfunction systems discussed by Reuter [123] and recurrent IFSes intro-duced recently by Barnsley [8]. The line segment J (F (2s)) is gener-ated using an IFS, for example consisting of two scaling transformationsQ1 and Q2 which map it onto its upper and lower half (Figure 8.7a).Subsequently, transformations T1, T2 and T3 are applied to create otherpoints of the set A(s). The order of transformation application is im-portant. Transformations Q1 and Q2 are used solely for the purposeof initial segment creation. They must not be applied after T1, T2 orT3, since in this case they would affect the branching structure underconsideration. The admissible sequences of transformations can be de-fined using a directed control graph (Figure 8.7b), and correspond tothe infinite set of paths starting at node a.1 The term controlled it-erated function system (CIFS) denotes an IFS with restrictions on thetransformation sequences imposed by a control graph. Thus, notingthe angle increment associated with symbols + and − by δ, the fractal

1Formally, the sequence of admissible transformations is the regular languageaccepted by the finite (Rabin-Scott) automaton represented by the graph in Fig-ure 8.7b.

8.2. Plant models and iterated function systems 185

Figure 8.7: Construction of the set A(S): (a) definition of an IFS {Q1, Q2}that generates the initial line segment, (b) the control graph specifying theadmissible sequences of transformation application

approximation of the leaf in Figure 5.11a is given by the CIFS with the Resulting CIFScontrol graph in Figure 8.7b and the transformations specified below:

Q1 =

0.5 0 0

0 0.5 00 0 1

Q2 =

0.5 0 0

0 0.5 00 s 1

T1 =

1/R cos δ 1/R sin δ 0

−1/R sin δ 1/R cos δ 00 s 1

T2 =

1/R cos δ −1/R sin δ 0

1/R sin δ 1/R cos δ 00 s 1

T3 =

1/R 0 0

0 1/R 00 2s 1

The second approach to the generation of the line segment J (F (2s)) Noninvertibletransforma-tions

is consistent with the method applied by Barnsley to specify the fernleaf in Figure 8.2. The idea is to map the entire branching structureA(s) onto the line J (F (2s)). This can be achieved using a noninvertibletransformation Q which collapses all branches into a vertical line. Thescaling factor along the y axis is the ratio of the desired segment length2s, and the limit height of the entire structure A(s),

h =2s

1 − 1/R.

186 Chapter 8. Fractal properties of plants

Figure 8.8: Two renderings of the compound leaf from Figure 5.11a, gener-ated using iterated function systems

This last value is calculated as the limit of the geometric series with thefirst term equal to 2s and the ratio equal to 1/R. Thus, the compoundleaf of Figure 5.11a is defined by an IFS consisting of transformation

Q =

0 0 0

0 1 − 1/R 00 0 1

and transformations T1, T2 and T3 specified as in the case of the con-trolled IFS.

Two fractal-based renderings of the set A(s) are shown in Fig-Renderingexamples ure 8.8. Figure 8.8a was obtained using the controlled IFS and a deter-

ministic algorithm to traverse the tree of admissible transformations.Figure 8.8b was obtained using the “ordinary” IFS and the randomselection of transformations. Figure 8.9 shows another fractal-basedrendering of the same structure. The spheres have radii equal to thedistance from the sphere center to the leaf, within a specified ε.

Figure 8.9: Carrot leaf �

8.2. Plant models and iterated function systems 187

188 Chapter 8. Fractal properties of plants

L-system with elongating internodes

‖L-system transformation

⇓L-system with decreasing apices

‖L-system analysis

⇓Recurrent equation in the domain of strings

‖Graphical interpretation

⇓Recurrent equation in the domain of sets

‖Passage to limit

⇓An equation expressing the limit set as a union ofthe limit object and reduced copies of itself

‖Equation solution

⇓An equation expressing the limit object as the im-age of an initial object under an iteration of aunion of transformations

‖Elimination of the initial object

⇓A (controlled) IFS

Figure 8.10: Steps in the construction of an IFS given an L-system capturinga developmental model

8.2. Plant models and iterated function systems 189

It is instructive to retrace the logical construction that started with Conclusionsan L-system, and ended with an iterated function system which cangenerate fractal approximations of the same object (Figure 8.10). Ananalysis of the operations performed in the subsequent steps of thisconstruction reveals its limitations, and clarifies the relationship be-tween strictly self-similar structures and real plants. The critical stepis the transformation of the L-system with elongating internodes to theL-system with decreasing apices. It can be performed as indicated inthe example if the plant maintains constant branching angles as well asfixed proportions between the mother and daughter segments, indepen-dent of branch order. This, in turn, can be achieved if all segments inthe modeled plant elongate exponentially over time. These are strongassumptions, and may be satisfied to different degrees in real plants.Strict self-similarity is an abstraction that captures the essential prop-erties of many plant structures and represents a useful point of referencewhen describing them in detail.

Epilogue

This quiet place, reminiscient of Claude Monet’s 1899 painting Water-lilies pool — Harmony in green, does not really exist. The scene wasmodeled using L-systems that captured the development of trees andwater plants, and illuminated by simulated sunlight. It is difficult notto appreciate how far the theory of L-systems and the entire field ofcomputer graphics have developed since their beginnings in the 1960’s,making such images possible. Yet the results contained in this book arenot conclusive and constitute only an introduction to the research onplant modeling for biological and graphics purposes. The algorithmicbeauty of plants is open to further exploration.

� Figure E.1: Water-lilies

Appendix A

Software environment forplant modeling

This book is illustrated with images of plants which exist only as math-ematical models visualized by means of computer graphics. The soft-ware environment used to construct and experiment with these modelsincludes dozens of programs and hundreds of data files. This createsthe nontrivial problem of organizing all components for easy definition,saving, retrieval and modification of the models. In order to solve it, theidea of simulation was extended beyond the level of individual plantsto an entire laboratory in botany [98]. Thus, a user can create andconduct experiments in a virtual laboratory by applying intuitive con-cepts and techniques from the “real” world. As an operating systemdefines the way a user perceives a computing environment, the virtuallaboratory determines a user’s perception of the environment in whichsimulated experiments take place. In the future, a virtual laboratorymay complement, extend, or even replace books as a means for gather-ing and presenting scientific information. Because of this potential, thelaboratory in which the research reported in this book was produced isdescribed here in more detail.

A.1 A virtual laboratory in botany

A virtual laboratory, like its “real” counterpart, is a playground for ex- User’sperspectiveperimentation. It comes with a set of objects pertinent to its scientific

domain (in this case, plant models), tools which operate on these ob-jects, a reference book and a notebook. Once the concepts and tools areunderstood, the user can expand the laboratory by adding new objects,creating new experiments, and recording descriptions in the notebook.An experienced user can expand the laboratory further by creating andinstalling new tools.

194 Appendix A. Software environment for plant modeling

Technically, a virtual laboratory is a microworld which can be ex-Laboratory =microworld +hypertext

plored under the guidance of a hypertext system. The term “microworld”denotes an interactive environment for creating and conducting simu-lated experiments. The guidance could be provided in the form of atraditional book, but an electronic document is more suitable for inte-gration with a microworld. In a sense, both components of the virtuallaboratory are described by Nelson in Dream Machines [104]. The pi-oneering role of this book in introducing the concept of hypertext isknown, but under the heading The Mind’s Eye the notion of a mi-croworld is also anticipated:

Suppose that you have a computer.What sorts of things would you do with it?Things that are imaginative

and don’t require too much else.I am hinting at something.You could have it make pictures and show you stuffand change what it shows depending on what you do.

A virtual laboratory can be divided into two components: the ap-Requirementsplication programs, data files and textual descriptions that comprisethe experiments; and the system support that provides the frameworkon which these domain-dependent experiments are built. The followinglist specifies the features of this framework.

• Consistent organization of the lab. In the lab environment,experiments are run by applying tools (programs) to objects (datafiles). An object consists of files that are grouped together sothat they can be retrieved easily. The format of the objects issufficiently standardized to allow straightforward implementationof common operations such as object saving and deletion.

• Inheritance of features. It is often the case that several ob-jects differ only in details. For example, two lilac inflorescencesmay differ only in the color of their petals. The mechanism ofinheritance is employed to store such objects efficiently.

• Version control. Interaction with an object during experimen-tation may result in a temporary or permanent modification. Inthe latter case, the user is able to decide whether the newly cre-ated object replaces the old one or should be stored as anotherversion of the original object.

• Interactive manipulation of objects. The laboratory pro-vides a set of general-purpose tools for manipulating object pa-rameters. For example, objects can be modified using controlpanels or by editing specific fields in a textual description of anexperiment.

A.1. A virtual laboratory in botany 195

• Flexibility in conducting experiments. The user may applytools to objects in a dynamic way while an experiment is beingconducted. This can be contrasted to a static experiment de-signed when the object is initially incorporated into the system.

• Guidance through the laboratory. A hypertext system im-poses a logical organization on the set of objects, provides atextual description of the experiments, and makes it possible tobrowse through the experiments in many ways. Specific experi-ments are invoked automatically when the corresponding text isselected, in order to facilitate demonstrations and assist a noviceuser.

So far, objects have been referred to in an intuitive way, relying Objectson the analogy between a real and virtual laboratory. For example, ifour interest is in the development of the gametophyte Microsorium lin-guaeforme, in a real laboratory we would experiment with a specimenof the plant, while in a virtual laboratory we explore the correspondingmathematical model. However, the analogy to real objects does not ex-tend to the level of detailed object definition. Specific design decisionsare needed for software development purposes. In the current design, alaboratory object is defined as a directory containing two types of filesand a subdirectory.

• The data files comprise our knowledge of a particular model.

• A specification file defines the data files which make up the objectand the tools which apply to them.

• A directory of extensions lists objects which inherit some featuresof the current object.

The object-oriented file structure which provides the basis for lab op-eration can be represented by a hierarchy of directories and files (Fig-ure A.1).

The path of subdirectories leading to an object establishes the in- Inheritance offeaturesheritance structure for the lab. Inheritance is based on the idea of

specifying new objects in reference to objects which already exist [81].The “old” object is called a prototype and the new one is its extension.The extension contains only those files which are different from thecorresponding files in the prototype. Files that remain the same aredelegated to the prototype by establishing links. In other words, theobject directory will contain those files that are unique to the object,and links to files that are inherited from its prototype (Figure A.2).This approach saves space, facilitates creation of objects similar to theprototype, and allows a single change in the prototype to propagatethrough all descendents.

196 Appendix A. Software environment for plant modeling

Figure A.1: The hierarchical structure of objects

Figure A.2: A prototype and its extension. Shaded areas indicate links.

A.1. A virtual laboratory in botany 197

Figure A.3: An object icon with menus

To conduct an experiment, all files that make up the selected object Version controlare copied to a temporary location called the lab table. Consequently,manipulation of object parameters does not disturb the stored version.When the experiment is finished, the user may save the results byoverwriting the original object or by creating an extension. In thelatter case, the files on the lab table are compared with those in theprototype object; those files that differ from the prototype are saved,and links to the remaining files are established automatically.

The ability to manipulate the parameters in an experiment easily is Objectmanipulationan essential feature of the virtual laboratory. As a rule, all parameters

involved in an experiment are supplied to the tools through the object’sdata files. In order to modify a parameter, the user edits the appropri-ate file, which is subsequently re-read by the application. Though theediting of parameters can be accomplished using a text editor, in manycases parameter modification can be performed more conveniently usingvirtual control panels [114]. The current implementation of the labora-tory provides the user with a general-purpose control manager whichcreates panels according to user-supplied configuration files.

The user is able to apply a tool to an object as a whole, without Toolapplicationdetailed knowledge of the programs involved or the component files.

This is achieved through the object’s specification file which lists allfiles associated with an object and the tools that can be applied tothem. This information is used to create a hierarchy of menus associ-ated with an icon representing the object (Figure A.3). The end nodesin the hierarchy invoke tools that operate on the object. For example,selection of the item image followed by the item generate from themenus in the figure would invoke the plant modeling program Pfg.

A user may browse through the objects in the lab by following either Browsingthe hierarchical structure of objects or hypertext links. The browser isused to navigate through the hierarchy, moving down through succes-sive extensions or up through previous levels. At any time, the usermay request that an object be placed on the lab table. The hypertextdocument associated with the lab provides an alternative method ofbrowsing and a means of relating objects independent of the hierarchy.

198 Appendix A. Software environment for plant modeling

A.2 List of laboratory programs

The essential programs incorporated into the virtual laboratory in botanyare listed below.

• Plant and fractal generator (Pfg)P. Prusinkiewicz and J. HananGiven an L-system, a set of viewing parameters and optional filesspecifying predefined surfaces, Pfg generates the modeled struc-ture by carrying out the derivation, then interpreting the resultingstring using turtle geometry. Both non-parametric and paramet-ric L-systems are supported. The model can be visualized directlyon the screen of an IRIS workstation or output to a file. The firstmode of operation is used to experiment with the model interac-tively and present developmental sequences. The output file canbe either in Postscript format, particularly suitable for printingresults such as fractal curves and inflorescence diagrams on a laserprinter, or in the format required by the ray-tracer Rayshade forrealistic rendering of the modeled structures.

• Modeling program for phyllotactic patterns (Spiral)D. R. FowlerSpiral is an interactive program for modeling organs with spiralphyllotactic patterns. The user can choose between planar andcylindrical patterns, and modify parameters which define modelgeometry (Chapter 4). This technique is faster than “growing”organs using parametric L-systems. Once an organ has been de-signed, it can be expressed using an L-system and incorporatedinto a plant structure.

• Interactive surface editor (Ise)J. HananIse makes it possible to define and modify bicubic surfaces con-sisting of one or several arbitrarily connected patches. The outputfiles produced by Ise are compatible with Pfg and Spiral.

• Modeling program for cellular structures (Mapl)F. D. FracchiaMapl accepts the specification of a two-dimensional cell layer cap-tured by a map L-system and generates the resulting developmen-tal sequence using the dynamic method of map interpretation.Options include map generation on the surface of a sphere, andthe simulation of development in three dimensions according toa given cellwork L-system. As in the case of Pfg, the models canbe visualized directly on the screen or output to a file in eitherPostscript or Rayshade format.

A.2. List of laboratory programs 199

• Control panel manager (Panel)L. Mercer and A. SniderThis program creates control panels containing sliders and but-tons, according to a configuration file provided by the user. Uponactivation of a control by the mouse, Panel generates a messagewhich indicates the corresponding control value. Application pro-grams process this information and modify the appropriate pa-rameters. For example, a panel can be used to control parametersused by Pfg, Spiral, or Mapl.

• Ray tracer (Rayshade)C. Kolb, Yale University

Rayshade reads a scene description from a text file, and renders itusing ray tracing. Scenes can be composed of primitives such asplanes, triangles, polygons, spheres, cylinders, cones and heightfields, grouped together to form objects. These objects can beinstantiated in other object definitions to create a hierarchicaldescription of a scene. Transformations including translation, ro-tation and scaling, and a variety of procedural textures can be ap-plied to any object. Extended light sources, simulation of depthof field, and adaptive supersampling are supported. The programuses 3D grids to partition object space for fast intersection tests.

• Previewer for the ray tracer (Preray)A. SniderPreray is a previewer for Rayshade used to provide a fast wireframe rendering of a scene before committing time to ray tracing.A control panel associated with Preray makes it possible to setviewing parameters interactively.

• Modeling program based on Euclidean constructions (L.E.G.O.)N. FullerL.E.G.O. makes it possible to model two- and three-dimensionalobjects using geometric constructions. In the scope of this book,L.E.G.O. was used to model man-made objects such as the Zinniavase and the Water-lilies bridge.

• Iterated function system generator (Ifsg)D. Hepting

A fractal defined by an iterated function system is described bya finite set of contractive affine transformations with an optionalfinite state control mechanism. Ifsg accepts input from a file spec-ifying the transformations and rendering information. The pro-gram is capable of rendering by either attracting, distance-basedor escape-time methods. The output can be displayed directlyon an IRIS workstation or written to a file for further processing.In the scope of this book, Ifsg was used to obtain results whichrelated plant models expressed using L-systems to fractals.

200 Appendix A. Software environment for plant modeling

Figure A.4: A virtual laboratory screen

Figure A.4 presents a sample screen of a Silicon Graphics IRIS 4D/60workstation running some of the above programs within the virtuallaboratory framework. The icon in the top right corner represents thelaboratory browser which was used to select a sunflower plant as thecurrent object. The icon underneath and the associated menu weresubsequently applied to select tools which operate on the object. Thecontrol panel in the bottom right corner of the screen is a part of thesurface editor Ise. The manipulated petal is displayed as a wire framein the window labeled Ise, and incorporated into a flower head by themodeling program Spiral which presents its output in the window sun-flower. The flower heads are in turn incorporated into a complete plantmodel generated by Pfg and rendered using Rayshade in the windowplant.rle. The panel below that window makes it possible to chooseorgans included in the model and change parameters related to theangles of the branching structure. The metaphor of a virtual labora-tory provides a uniform interface to various operations on the selectedplant, ranging from the modification of a petal to the rendering of thecomplete model.

Appendix B

About the figures

The following descriptions of the color images include details about thepictures not described in the main text. Unless otherwise stated, figureswere created at the University of Regina.

Figure 1.19 [page 20] Three-dimensional Hilbert curveF. D. Fracchia, P. Prusinkiewicz, N. Fuller(1989)

This image was rendered using ray-tracing without shadows.

Figure 1.25 [page 26] Three-dimensional bushP. Prusinkiewicz (1986)

Simple branching structure, rendered using the firmware of a Sil-icon Graphics IRIS workstation. Total generating and renderingtime on IRIS 4D/20: 4 seconds.

Figure 1.28 [page 29] Flower fieldP. Prusinkiewicz (1986)

The field contains four rows of four plants. The scene was ren-dered with IRIS firmware, using depth-cueing to assign colors topetals.

Figure 1.35 [page 45] Developmental stages of Anabaena catenulaJ. Hanan, P. Prusinkiewicz (1989)

Figure 2.1 [page 52] Organic architectureNed Greene, NYIT (1989)

202 Appendix B. About the figures

An array of 300 x 300 x 300 voxel space automata was used totrack a polygonal model of a house. Rendering was performedusing a probabilistic radiosity method. See [54] for a full descrip-tion.

Figure 2.3 [page 54] Acer graphicsJules Bloomenthal, NYIT (1984)

A model of a maple tree. The basic branching structure was gen-erated recursively. Limbs were modeled as generalized cylinders,obtained by moving discs of varying radii along spline curves.Real bark texture was digitized and used as a bump map. Leaftexture was obtained by digitizing a photograph of a real leaf andemphasizing the veins using a paint program. See [11] for details.

Figure 2.4 [page 54] Forest sceneBill Reeves, Pixar (1984)

A scene from the film The Adventures of Andre and Wally B,modeled using particle systems. Shading and shadows were ap-proximated using probabilistic techniques. Visible surfaces weredetermined using depth-sorting. See [119] for a full description.

Figure 2.5 [page 55] Oil palm tree canopyCIRAD Modelisation Laboratory (1990)

A developmental model of oil palm trees, modeled using themethod originated by de Reffye and described from the graph-ics perspective in [30].

Figure 2.10 [page 61] Medicine lakeF. K. Musgrave, C. E. Kolb, P. Prusinkiewicz,B. B. Mandelbrot (1988)

A scene combining a fractal terrain model, a tree generated us-ing L-systems, and a rainbow. The rainbow model was derivedfrom a simulation of refraction with dispersion of light throughan idealized raindrop. Procedural textures were applied to themountains, the water surface and a vertical plane modeling thesky. See [101] for further details.

Figure 2.11 [page 62] Surrealistic elevatorA. Snider, P. Prusinkiewicz, N. Fuller (1989)

The elevator was modeled using L.E.G.O. The island is a su-perquadratic surface. Procedural textures were applied to createstars in the sky, craters on the moon, colored layers in the rock,

Appendix B. About the figures 203

waves in the lake and imperfections in the glass that covers theelevator.

Figure 3.2 [page 69] CrocusesJ. Hanan, D. R. Fowler (1990)

The petals were modeled as Bezier surfaces, with the shapes de-termined using Ise.

Figure 3.4 [page 72] Lily-of-the-valleyP. Prusinkiewicz, J. Hanan (1987)

Figure 3.5 [page 74] Development of Capsella bursa-pastorisP. Prusinkiewicz, A. Lindenmayer (1987)

Figure 3.6 [page 75] Apple twigP. Prusinkiewicz, D. R. Fowler (1990)

This twig model was developed in one spring day, looking at areal twig nearby. This time is indicative for most inflorescencemodels shown.

Figure 3.11 [page 81] A mintP. Prusinkiewicz (1988)

Figure 3.14 [page 84] Development of Lychnis coronariaP. Prusinkiewicz, J. Hanan (1987)

Figure 3.17 [page 90] Development of Mycelis muralisP. Prusinkiewicz, A. Lindenmayer (1987)

Figure 3.18 [page 91] A three-dimensional rendering of the MycelismodelsP. Prusinkiewicz, J. Hanan (1987)

All internodes in the model are assumed to have the same length.In reality, the internodes have different lengths, and the structureis less crowded.

Figure 3.19 [page 92] Lilac inflorescencesP. Prusinkiewicz, J. Hanan, D. R. Fowler(1990)

204 Appendix B. About the figures

Figure 3.21 [page 94] The Garden of LP. Prusinkiewicz, F. D. Fracchia, J. Hanan,D. R. Fowler (1988)

All plants were modeled with L-systems and rendered using theIRIS firmware. Images corresponding to different viewing planes(the background lilac twigs, the apple twig and the daisies) weredefocused separately using low-pass filters to simulate the depthof field, then composited with a focused image of lilac inflores-cences. The sky was generated using a fractal algorithm.

Figure 3.23 [page 96] Wild carrotP. Prusinkiewicz (1988)

Figure 4.3 [page 102] Close-up of a daisy capitulumD. R. Fowler (1988)

The petals and florets were modeled as Bezier surfaces.

Figure 4.4 [page 102] Domestic sunflower headD. R. Fowler, P. Prusinkiewicz (1989)

Figure 4.5 [page 105] Sunflower fieldD. R. Fowler, N. Fuller, J. Hanan, A. Snider(1990)

This image contains approximately 400 plants, each with 15 flow-ers. A flower has 21 petals and 300 seeds, modeled using 600 tri-angles and 400 triangles respectively. Counting leaves and buds,the entire scene contains about 800,000,000 triangles. The imagewas ray-traced with adaptive supersampling on a grid of 1024 x768 pixels using 45 hours of CPU time on a MIPS M-120 com-puter.

Figure 4.6 [page 106] ZinniasD. R. Fowler, P. Prusinkiewicz, J. Hanan,N. Fuller (1990)

The vase was modeled using L.E.G.O. and rendered with a pro-cedural texture. The scene was illuminated by one extended lightsource.

Appendix B. About the figures 205

Figure 4.7 [page 106] Close-up of zinniasD. R. Fowler, P. Prusinkiewicz, A. Snider(1990)

This scene was rendered using distributed ray-tracing to simu-late the depth field.

Figure 4.8 [page 108] Water-lilyD. R. Fowler, J. Hanan (1990)

Figure 4.9 [page 108] Lily pondD. R. Fowler, J. Hanan, P. Prusinkiewicz,N. Fuller (1990)

The wavelets on the water surface were obtained using bump-mapping with a procedurally defined texture.

Figure 4.10 [page 109] RosesD. R. Fowler, J. Hanan, P. Prusinkiewicz(1990)

Distributed ray-tracing with one extended light source was usedto simulate depth of field and create fuzzy shadows.

Figure 4.11 [page 111] Parastichies on a cylinderD. R. Fowler (1990)

Figure 4.15 [page 116] PineapplesD. R. Fowler, A. Snider (1990)

The image incorporates a physically-based model of a tableclothapproximated as an array of masses connected by springs andplaced in a gravitational field. The scene is illuminated by threeextended light sources.

Figure 4.16 [page 117] Spruce conesD. R. Fowler, J. Hanan (1990)

Figure 4.17 [page 117] Carex laevigataJ. Hanan, P. Prusinkiewicz (1989)

The entire plant, including the leaves, was modeled using para-metric L-systems.

206 Appendix B. About the figures

Figure 5.2 [page 121] Maraldi figureNed Greene, NYIT (1984)

The shapes of leaves, calyxes and petals were defined using a paintprogram, by interpreting gray levels as height. Painted textureswere mapped onto the surfaces of leaves and calyxes. Smoothgradation of color across the petals was obtained by assigningcolors to the vertices of the polygon meshes representing flowers,then interpolating colors across polygons using Gouraud shading.The vines were rendered with bump-mapping, using a digitizedimage of real bark.

Figure 5.3 [page 121] The fernP. Prusinkiewicz (1986)

Figure 5.7 [page 125] A rose in a vaseD. R. Fowler, J. Hanan, P. Prusinkiewicz(1990)

Petals and thorns are Bezier surfaces incorporated into a rosemodel expressed using L-systems. The vase was modeled as asurface of revolution.

Figure 6.3 [page 141] Development of Anabaena catenulaP. Prusinkiewicz, F. D. Fracchia (1989)

Each developmental stage is plotted in one scan line.

Figure 7.13 [page 161] Simulated development of MicrosoriumlinguaeformeF. D. Fracchia, P. Prusinkiewicz,M. J. M. de Boer (1989)

Cells are represented as polygons, rendered using the IRIS firmware.The development can be visualized directly on the screen of anIRIS 4D/20 workstation without resorting to single-frame anima-tion techniques.

Figure 7.14 [page 161] Microphotograph of Microsorium linguaeformeM. J. M. de Boer, University of Utrecht

Figure 7.16 [page 163] Simulated development of DryopteristhelypterisF. D. Fracchia, P. Prusinkiewicz,M. J. M. de Boer (1989)

Appendix B. About the figures 207

Figure 7.19 [page 169] Developmental sequence of Patella vulgataF. D. Fracchia, A. Lindenmayer,M. J. M. de Boer (1989)

Cells are represented as spheres. Intersections of spheres insidethe modeled embryo are ignored, since they do not affect theray-traced images.

Figure 7.20 [page 169] An electron microscope image of PatellavulgataW. J. Dictus, University of Utrecht

Figure 8.4 [page 180] Fern duneP. Prusinkiewicz, D. Hepting (1989)

The shape of the leaf has been captured using a controlled iter-ated function system. A continuous escape-time function definespoint altitudes, resulting in a surrealistic incorporation of a leafinto the landscape.

Figure 8.9 [page 186] Carrot leafD. Hepting, P. Prusinkiewicz (1989)

The leaf shape has been modeled using a controlled iterated func-tion system. The scene consists of a set of spheres, with the radiusequal to the distance to the leaf. The image was rendered usingray-tracing.

Figure E.1 [page 191] Water-liliesD. R. Fowler, J. Hanan, P. Prusinkiewicz,N. Fuller (1990)

A scene inspired by Water-lilies pool - Harmony in green byClaude Monet (1899). All trees and water-lilies were modeled us-ing L-systems. The willow twigs bend downwards due to a strongtropism effect, simulating gravity. The bridge was modeled usingL.E.G.O. The sky is a sphere with a procedural texture. Theentire scene was ray-traced, then the resulting image was repre-sented as a set of small circles, with the colors close but not equalto the average of pixel colors underneath. This last operation wasaimed at creating the appearance of an impressionistic painting.

Figure A.4 [page 200] Virtual labL. Mercer, D. R. Fowler (1990)

Turtle interpretation ofsymbols

Symbol Interpretation Page

F Move forward and draw a line. 7, 46

f Move forward without drawing a line. 7, 46

+ Turn left. 7, 19, 46

− Turn right. 7, 19

∧ Pitch up. 19, 46

& Pitch down. 19, 46

\ Roll left. 19, 46

/ Roll right. 19, 46

| Turn around. 19, 46

$ Rotate the turtle to vertical. 57

[ Start a branch. 24

] Complete a branch. 24

{ Start a polygon. 120, 127

G Move forward and draw a line. Do not record a vertex. 122

. Record a vertex in the current polygon. 122, 127

} Complete a polygon. 120, 127

∼ Incorporate a predefined surface. 119

! Decrement the diameter of segments. 26, 57′ Increment the current color index. 26

% Cut off the remainder of the branch. 73

Bibliography

[1] H. Abelson and A. A. diSessa. Turtle geometry. M.I.T. Press,Cambridge, 1982.

[2] M. Aono and T. L. Kunii. Botanical tree image generation. IEEEComputer Graphics and Applications, 4(5):10–34, 1984.

[3] M. J. Apter. Cybernetics and development. Pergamon Press,Oxford, 1966. (International Series of Monographs in Pure andApplied Biology/Zoology Division Vol. 29).

[4] W. W. Armstrong. The dynamics of tree linkages with a fixedroot link and limited range of rotation. Actes du Colloque Inter-nationale l’Imaginaire Numerique ’86, pages 16–21, 1986.

[5] J. W. Backus. The syntax and semantics of the proposed interna-tional algebraic language of the Zurich ACM-GAMM conference.In Proc. Intl. Conf. on Information Processing, pages 125–132.UNESCO, 1959.

[6] B. I. Balinsky. An introduction to embryology. W. B. Saunders,Philadelphia, 1970.

[7] M. F. Barnsley. Fractals everywhere. Academic Press, San Diego,1988.

[8] M. F. Barnsley, J. H. Elton, and D. P. Hardin. Recurrent iteratedfunction systems. Constructive Approximation, 5:3–31, 1989.

[9] B. A. Barsky. The Beta-spline: A local representation based onshape parameters and fundamental geometric measures. PhD the-sis, Department of Computer Science, University of Utah, 1981.

[10] R. Bartels, J. Beatty, and B. Barsky, editors. An introductionto splines for use in computer graphics and geometric modeling.Morgan Kaufman, Los Altos, California, 1987.

[11] J. Bloomenthal. Modeling the mighty maple. Proceedings ofSIGGRAPH ’85 (San Francisco, California, July 22-26, 1985) inComputer Graphics, 19, 3 (July 1985), pages 305–311, ACM SIG-GRAPH, New York, 1985.

212 Bibliography

[12] B. G. Briggs and L. A. S. Johnson. Evolution in the Myrtaceae –evidence from inflorescence structure, Appendix I: The relevanceof Troll’s system of inflorescence typology. In Proceedings of theLinnean Society of New South Wales, volume 102, pages 236–240,1979.

[13] N. Chomsky. Three models for the description of language. IRETrans. on Information Theory, 2(3):113–124, 1956.

[14] V. Claus, H. Ehrig, and G. Rozenberg, editors. Graph grammarsand their application to computer science; First InternationalWorkshop. Lecture Notes in Computer Science 73. Springer-Verlag, Berlin, 1979.

[15] D. Cohen. Computer simulation of biological pattern generationprocesses. Nature, 216:246–248, 1967.

[16] E. Costes. Analyse architecturale et modelisation du litchi. PhDthesis, Universite des Sciences et Techniques du Languedoc, 1988.

[17] H. S. M. Coxeter. Introduction to geometry. J. Wiley & Sons,New York, 1961.

[18] H. S. M. Coxeter. The role of intermediate convergents in Tait’sexplanation for phyllotaxis. J. Algebra, 20:167–175, 1972.

[19] K. Culik II and D. Wood. A mathematical investigation of prop-agating graph OL-systems. Information and Control, 43:50–82,1979.

[20] P. Dabadie, P. de Reffye, and P. Dinouard. Modelisation de lacroissance et de l’architecture d’un bambou. In Deuxieme congresinternational du bambou, 1988.

[21] C. Davis and D. E. Knuth. Number representations and dragoncurves. J. of Recreational Mathematics, 3:66–81, 133–149, 1970.

[22] M. J. M. de Boer. Analysis and computer generation of divi-sion patterns in cell layers using developmental algorithms. PhDthesis, University of Utrecht, the Netherlands, 1989.

[23] M. J. M. de Boer and A. Lindenmayer. Map 0L-systems withedge label control: Comparison of marker and cyclic systems.In H. Ehrig, M. Nagl, A. Rosenfeld, and G. Rozenberg, editors,Graph-grammars and their application to computer science, Lec-ture Notes in Computer Science 291, pages 378–392. Springer-Verlag, 1987.

Bibliography 213

[24] M. de Does and A. Lindenmayer. Algorithms for the genera-tion and drawing of maps representing cell clones. In H. Ehrig,M. Nagl, and G. Rozenberg, editors, Graph-grammars and theirapplication to computer science, Lecture Notes in Computer Sci-ence 153, pages 39–57. Springer-Verlag, 1983.

[25] C. G. de Koster and A. Lindenmayer. Discrete and continuousmodels for heterocyst differentiation in growing filaments of blue-green bacteria. In Acta Biotheoretica, volume 36, pages 249–273.Kluwer Academic, the Netherlands, 1987.

[26] P. de Reffye. Modele mathematique aleatorie et simulation dela croissance et de l’architecture du cafeier robusta. Premierepartie, Cafe-Cacao-The, 25(2):83–104, 1981. Deuxieme partie,Cafe-Cacao-The, 25(4):219–230, 1981. Troisieme partie, Cafe-Cacao-The, 26(2):77–96, 1982. Quatrieme partie, Cafe-Cacao-The, 27(1):3–20, 1983.

[27] P. de Reffye. Travaux sur la cafeier. CIRAD, Montpellier. Col-lection of earlier papers.

[28] P. de Reffye. Modelisation et simulation de la verse de cafeier, al’aide de la theorie de la resistance des materiaux. Cafe-Cacao-The, XX(4):251–272, 1976.

[29] P. de Reffye, M. Cognee, M. Jaeger, and B. Traore. Modelisationde la croissance et de l’architecture du cotonnier. Manuscript,1988.

[30] P. de Reffye, C. Edilin, J. Francon, M. Jaeger, and C. Puech.Plant models faithful to botanical structure and development.Proceedings of SIGGRAPH ’88 (Atlanta, Georgia, August 1-5,1988), in Computer Graphics 22,4 (August 1988), pages 151–158,ACM SIGGRAPH, New York, 1988.

[31] F. M. Dekking. Recurrent sets. Advances in Mathematics,44(1):78–104, 1982.

[32] F. M. Dekking. Recurrent sets: A fractal formalism. Report82-32, Delft University of Technology, 1982.

[33] H. Ehrig, M. Nagl, A. Rosenfeld, and G. Rozenberg, editors.Graph grammars and their application to computer science; ThirdInternational Workshop. Lecture Notes in Computer Science 291.Springer-Verlag, Berlin, 1987.

[34] H. Ehrig, M. Nagl, and G. Rozenberg, editors. Graph grammarsand their application to computer science; Second InternationalWorkshop. Lecture Notes in Computer Science 153. Springer-Verlag, Berlin, 1983.

214 Bibliography

[35] P. Eichhorst and W. J. Savitch. Growth functions of stochas-tic Lindenmayer systems. Information and Control, 45:217–228,1980.

[36] R. O. Erickson. The geometry of phyllotaxis. In J. E. Dale andF. L. Milthrope, editors, The growth and functioning of leaves,pages 53–88. University Press, Cambridge, 1983.

[37] G. Eyrolles. Synthese d’images figuratives d’arbres par desmethodes combinatoires. PhD thesis, Universite de BordeauxI, 1986.

[38] J. B. Fisher and H. Honda. Computer simulation of branchingpattern and geometry in Terminalia (Combretaceae), a tropicaltree. Botanical Gazette, 138(4):377–384, 1977.

[39] J. B. Fisher and H. Honda. Branch geometry and effective leafarea: A study of Terminalia–branching pattern, Parts I and II.American Journal of Botany, 66:633–655, 1979.

[40] J. D. Foley and A. Van Dam. Fundamentals of interactive com-puter graphics. Addison-Wesley, Reading, Massachusetts, 1982.

[41] L. Fox and D. F. Mayers. Numerical solution of ordinary differ-ential equations. Chapman and Hall, London, 1987.

[42] F. D. Fracchia, P. Prusinkiewicz, and M. J. M. de Boer. Visualiza-tion of the development of multicellular structures. In Proceedingsof Graphics Interface ’90, pages 267–277, 1990.

[43] H. Freeman. On encoding arbitrary geometric configurations.IRE Trans. Electronic. Computers, 10:260–268, 1961.

[44] D. Frijters. Mechanisms of developmental integration of Asternovae-angliae L. and Hieracium murorum L. Annals of Botany,42:561–575, 1978.

[45] D. Frijters. Principles of simulation of inflorescence development.Annals of Botany, 42:549–560, 1978.

[46] D. Frijters and A. Lindenmayer. A model for the growth andflowering of Aster novae-angliae on the basis of table (1,0)L-systems. In G. Rozenberg and A. Salomaa, editors, L Systems,Lecture Notes in Computer Science 15, pages 24–52. Springer-Verlag, Berlin, 1974.

[47] D. Frijters and A. Lindenmayer. Developmental descriptions ofbranching patterns with paracladial relationships. In A. Linden-mayer and G. Rozenberg, editors, Automata, languages, develop-ment, pages 57–73. North-Holland, Amsterdam, 1976.

Bibliography 215

[48] M. Gardner. Mathematical games: An array of problems that canbe solved with elementary mathematical techniques. ScientificAmerican, 216, 1967. 3:124–129 (March), 4:116–123 (April).

[49] M. Gardner. Mathematical games: The fantastic combinationsof John Conway’s new solitaire game “life”. Scientific American,223(4):120–123, October 1970.

[50] M. Gardner. Mathematical games: On cellular automata, self-reproduction, the Garden of Eden and the game “life”. ScientificAmerican, 224(2):112–117, February 1971.

[51] M. Gardner. Mathematical games – in which “monster” curvesforce redefinition of the word “curve”. Scientific American,235(6):124–134, December 1976.

[52] S. Ginsburg and H. G. Rice. Two families of languages related toALGOL. J. ACM, 9(3):350–371, 1962.

[53] H. Gravelius. Flusskunde. Goschen, Berlin, 1914.

[54] N. Greene. Voxel space automata: Modeling with stochasticgrowth processes in voxel space. Proceedings of SIGGRAPH ’89(Boston, Mass., July 31-August 4, 1989), in Computer Graph-ics 23,4 (August 1989), pages 175–184, ACM SIGGRAPH, NewYork, 1989.

[55] B. E. S. Gunning. Microtubules and cytomorphogenesis in adeveloping organ: The root primordium of Azolla pinnata. InO. Kiermayer, editor, Cytomorphogenesis in plants, Cell BiologyMonographs 8, pages 301–325. Springer-Verlag, Wien, 1981.

[56] A. Habel and H.-J. Kreowski. On context-free graph languagesgenerated by edge replacement. In H. Ehrig, M. Nagl, andG. Rozenberg, editors, Graph grammars and their application tocomputer science; Second International Workshop, Lecture Notesin Computer Science 153, pages 143–158. Springer-Verlag, Berlin,1983.

[57] A. Habel and H.-J. Kreowski. May we introduce to you: Hy-peredge replacement. In H. Ehrig, M. Nagl, G. Rozenberg, andA. Posenfeld, editors, Graph grammars and their application tocomputer science; Third International Workshop, Lecture Notesin Computer Science 291, pages 15–26. Springer-Verlag, Berlin,1987.

[58] F. Halle, R. A. A. Oldeman, and P. B. Tomlinson. Tropical treesand forests: An architectural analysis. Springer-Verlag, Berlin,1978.

216 Bibliography

[59] P. H. Hellendoorn and A. Lindenmayer. Phyllotaxis in Bryophyl-lum tubiflorum: Morphogenetic studies and computer simula-tions. Acta Biol. Neerl, 23(4):473–492, 1974.

[60] D. Hepting, P. Prusinkiewicz, and D. Saupe. Rendering methodsfor iterated function systems. Manuscript, 1990.

[61] G. Herman, A. Lindenmayer, and G. Rozenberg. Description ofdevelopmental languages using recurrence systems. MathematicalSystems Theory, 8:316–341, 1975.

[62] G. T. Herman and G. Rozenberg. Developmental systems andlanguages. North-Holland, Amsterdam, 1975.

[63] D. Hilbert. Ueber stetige Abbildung einer Linie auf einFlachenstuck. Mathematische Annalin., 38:459–460, 1891.

[64] P. Hogeweg and B. Hesper. A model study on biomorphologicaldescription. Pattern Recognition, 6:165–179, 1974.

[65] H. Honda. Description of the form of trees by the parameters ofthe tree-like body: Effects of the branching angle and the branchlength on the shape of the tree-like body. Journal of TheoreticalBiology, 31:331–338, 1971.

[66] H. Honda and J. B. Fisher. Tree branch angle: Maximizing ef-fective leaf area. Science, 199:888–890, 1978.

[67] H. Honda and J. B. Fisher. Ratio of tree branch lengths: Theequitable distribution of leaf clusters on branches. Proceedings ofthe National Academy of Sciences USA, 76(8):3875–3879, 1979.

[68] H. Honda, P. B. Tomlinson, and J. B. Fisher. Computer simula-tion of branch interaction and regulation by unequal flow rates inbotanical trees. American Journal of Botany, 68:569–585, 1981.

[69] H. Honda, P. B. Tomlinson, and J. B. Fisher. Two geometricalmodels of branching of botanical trees. Annals of Botany, 49:1–11, 1982.

[70] R. E. Horton. Erosioned development of systems and theirdrainage basins, hydrophysical approach to quantitative morphol-ogy. Bull. Geol. Soc. America, 56:275–370, 1945.

[71] R. E. Horton. Hypsometric (area-altitude) analysis of erosionaltopology. Bull. Geol. Soc. America, 63:1117–1142, 1952.

[72] R. Hunt. Plant growth analysis. Studies in Biology 96. EdwardArnold, London, 1978.

Bibliography 217

[73] R. Hunt. Plant growth curves – the functional approach to plantgrowth analysis. Edward Arnold, London, 1982.

[74] J. E. Hutchinson. Fractals and self-similarity. Indiana UniversityJournal of Mathematics, 30(5):713–747, 1981.

[75] G. Van Iterson. Mathematische und mikroskopish-anatomischeStudien uber Blattstellungen. Gustav Fischer, Jena, 1907.

[76] M. Jaeger. Representation et simulation de croissance desvegetaux. PhD thesis, Universite Louis Pasteur de Strasbourg,1987.

[77] J. M. Janssen and A. Lindenmayer. Models for the control ofbranch positions and flowering sequences of capitula in Mycelismuralis (L.) Dumont (Compositae). New Phytologist, 105:191–220, 1987.

[78] R. V. Jean. Mathematical modelling in phyllotaxis: The state ofthe art. Mathematical Biosciences, 64:1–27, 1983.

[79] H. Jurgensen and A. Lindenmayer. Modelling development byOL-systems: Inference algorithms for developmental systemswith cell lineages. Bulletin of Mathematical Biology, 49(1):93–123, 1987.

[80] A. N. Kolmogorov. Three approaches to the quantitative defini-tion of information. Int. J. Comp. Math, 2:157–168, 1968.

[81] H. Lieberman. Using prototypical objects to implement sharedbehavior in object oriented systems. In Proceedings of theACM Conference on Object-Oriented Programming Systems,Languages, and Applications, pages 214–223, New York, 1986.Association for Computing Machinery.

[82] A. Lindenmayer. Mathematical models for cellular interactionin development, Parts I and II. Journal of Theoretical Biology,18:280–315, 1968.

[83] A. Lindenmayer. Adding continuous components to L-systems. InG. Rozenberg and A. Salomaa, editors, L Systems, Lecture Notesin Computer Science 15, pages 53–68. Springer-Verlag, Berlin,1974.

[84] A. Lindenmayer. Developmental algorithms: Lineage versus in-teractive control mechanisms. In S. Subtelny and P. B. Green,editors, Developmental order: Its origin and regulation, pages219–245. Alan R. Liss, New York, 1982.

218 Bibliography

[85] A. Lindenmayer. Models for plant tissue development with celldivision orientation regulated by preprophase bands of micro-tubules. Differentiation, 26:1–10, 1984.

[86] A. Lindenmayer. Positional and temporal control mechanisms ininflorescence development. In P. W. Barlow and D. J. Carr, ed-itors, Positional controls in plant development. University Press,Cambridge, 1984.

[87] A. Lindenmayer. An introduction to parallel map generating sys-tems. In H. Ehrig, M. Nagl, A. Rosenfeld, and G. Rozenberg,editors, Graph grammars and their application to computer sci-ence; Third International Workshop, Lecture Notes in ComputerScience 291, pages 27–40. Springer-Verlag, Berlin, 1987.

[88] A. Lindenmayer. Models for multicellular development: Char-acterization, inference and complexity of L-systems. In A. Kel-menova and J. Kelmen, editors, Trends, techniques and problemsin theoretical computer science, Lecture Notes in Computer Sci-ence 281, pages 138–168. Springer-Verlag, Berlin, 1987.

[89] A. Lindenmayer and P. Prusinkiewicz. Developmental modelsof multicellular organisms: A computer graphics perspective. InC. Langton, editor, Artificial Life: Proceedings of an Interdisci-plinary Workshop on the Synthesis and Simulation of Living Sys-tems held September, 1987, in Los Alamos, New Mexico, pages221–249. Addison-Wesley, Redwood City, 1989.

[90] A. Lindenmayer and G. Rozenberg, editors. Automata, languages,development. North-Holland, Amsterdam, 1976.

[91] A. Lindenmayer and G. Rozenberg. Parallel generation of maps:Developmental systems for cell layers. In V. Claus, H. Ehrig, andG. Rozenberg, editors, Graph grammars and their application tocomputer science; First International Workshop, Lecture Notesin Computer Science 73, pages 301–316. Springer-Verlag, Berlin,1979.

[92] J. Luck, A. Lindenmayer, and H. B. Luck. Models for celltetrads and clones in meristematic cell layers. Botanical Gazette,149:1127–141, 1988.

[93] J. Luck and H. B. Luck. Generation of 3-dimensional plant bod-ies by double wall map and stereomap systems. In H. Ehrig,M. Nagl, and G. Rozenberg, editors, Graph Grammars and TheirApplication to Computer Science; Second International Work-shop, Lecture Notes in Computer Science 153, pages 219–231.Springer-Verlag, Berlin, 1983.

Bibliography 219

[94] N. Macdonald. Trees and networks in biological models. J. Wiley& Sons, New York, 1983.

[95] B. B. Mandelbrot. The fractal geometry of nature. W. H. Free-man, San Francisco, 1982.

[96] D. M. McKenna. SquaRecurves, E-tours, eddies and frenzies:Basic families of Peano curves on the square grid. In Proceed-ings of the Eugene Strens Memorial Conference on RecreationalMathematics and its History, 1989. To appear.

[97] H. Meinhardt. Models of biological pattern formation. AcademicPress, New York, 1982.

[98] L. Mercer, P. Prusinkiewicz, and J. Hanan. The concept and de-sign of a virtual laboratory. In Proceedings of Graphics Interface’90, pages 149–155. CIPS, 1990.

[99] G. J. Mitchison and Michael Wilcox. Rules governing cell divisionin Anabaena. Nature, 239:110–111, 1972.

[100] D. Muller-Doblies and U. Muller-Doblies. Cautious improvementof a descriptive terminology of inflorescences. Monocot Newsletter4, 1987.

[101] F. K. Musgrave, C. E. Kolb, and R. S. Mace. The synthesisand rendering of eroded fractal terrains. Proceedings of SIG-GRAPH ’89 (Boston, Mass., July 31-August 4, 1989), in Com-puter Graphics 23,4 (August 1989), pages 41–50, ACM SIG-GRAPH, New York, 1989.

[102] A. Nakamura, A. Lindenmayer, and K. Aizawa. Some systemsfor map generation. In G. Rozenberg and A. Salomaa, editors,The Book of L, pages 323–332. Springer-Verlag, Berlin, 1986.

[103] P. Naur et al. Report on the algorithmic language ALGOL 60.Communications of the ACM, 3(5):299–314, 1960. Revised inComm. ACM 6(1):1-17.

[104] T. Nelson. Computer lib and dream machines. Self-published,1980.

[105] P. Oppenheimer. Real time design and animation of fractal plantsand trees. Computer Graphics, 20(4):55–64, 1986.

[106] G. Peano. Sur une courbe, qui remplit tout une aire plaine. Math.Annln., 36:157–160, 1890. Translated in G. Peano, Selected worksof Giuseppe Peano, H. C. Kennedy, editor, pages 143–149, Uni-versity of Toronto Press, Toronto, 1973.

220 Bibliography

[107] H. Peitgen and D. Saupe, editors. The science of fractal images.Springer-Verlag, New York, 1988.

[108] F. P. Preparata and R. T. Yeh. Introduction to Discrete Struc-tures. Addison-Wesley, Reading, Massachusetts, 1973.

[109] P. Prusinkiewicz. Graphical applications of L-systems. In Pro-ceedings of Graphics Interface ’86 — Vision Interface ’86, pages247–253. CIPS, 1986.

[110] P. Prusinkiewicz. Score generation with L-systems. In Proceedingsof the International Computer Music Conference ’86, pages 455–457, 1986.

[111] P. Prusinkiewicz. Applications of L-systems to computer imagery.In H. Ehrig, M. Nagl, A. Rosenfeld, and G. Rozenberg, editors,Graph grammars and their application to computer science; ThirdInternational Workshop, pages 534–548. Springer-Verlag, Berlin,1987. Lecture Notes in Computer Science 291.

[112] P. Prusinkiewicz and J. Hanan. Lindenmayer systems, frac-tals, and plants, volume 79 of Lecture Notes in Biomathematics.Springer-Verlag, Berlin, 1989.

[113] P. Prusinkiewicz and J. Hanan. Visualization of botanical struc-tures and processes using parametric L-systems. In D. Thalmann,editor, Scientific Visualization and Graphics Simulation, pages183–201. J. Wiley & Sons, 1990.

[114] P. Prusinkiewicz and K. Krithivasan. Algorithmic generation ofSouth Indian folk art patterns. In Proceedings of the Interna-tional Conference on Computer Graphics ICONCG ’88, Singa-pore, 1988.

[115] P. Prusinkiewicz, K. Krithivasan, and M. G. Vijayanarayana. Ap-plication of L-systems to algorithmic generation of South Indianfolk art patterns and karnatic music. In R. Narasimhan, editor,A perspective in theoretical computer science — commemorativevolume for Gift Siromoney, pages 229–247. World Scientific, Sin-gapore, 1989. Series in Computer Science Vol. 16.

[116] P. Prusinkiewicz, A. Lindenmayer, and F. D. Fracchia. Synthesisof space-filling curves on the square grid. To appear in Proceedingsof FRACTAL ’90, the 1st IFIP conference on fractals, Lisbon,Portugal, June 6-8, 1990.

[117] P. Prusinkiewicz, A. Lindenmayer, and J. Hanan. Developmen-tal models of herbaceous plants for computer imagery purposes.Proceedings of SIGGRAPH ’88 (Atlanta, Georgia, August 1-5,1988), in Computer Graphics 22,4 (August 1988), pages 141–150,ACM SIGGRAPH, New York, 1988.

Bibliography 221

[118] P. Prusinkiewicz and G. Sandness. Koch curves as attractors andrepellers. IEEE Computer Graphics and Applications, 8(6):26–40,1988.

[119] W. T. Reeves and R. Blau. Approximate and probabilistic al-gorithms for shading and rendering structured particle systems.Proceedings of SIGGRAPH ’85 (San Francisco, California, July22-26, 1985) in Computer Graphics, 19, 3 (July 1985), pages 313–322, ACM SIGGRAPH, New York, 1985.

[120] W. R. Remphrey, B. R. Neal, and T. A. Steeves. The morphologyand growth of Arctostaphylos uva-ursi (bearberry), parts i and ii.Canadian Journal of Botany, 61(9):2430–2458, 1983.

[121] W. R. Remphrey and G. R. Powell. Crown architecture of Larixlaricina saplings: Quantitative analysis and modelling of (non-sylleptic) order 1 branching in relation to development of themain stem. Canadian Journal of Botany, 62(9):1904–1915, 1984.

[122] W. R. Remphrey and G. R. Powell. Crown architecture of Larixlaricina saplings: Sylleptic branching on the main stem. Cana-dian Journal of Botany, 63(7):1296–1302, 1985.

[123] L. H. Reuter. Rendering and magnification of fractals using in-terated function systems. PhD thesis, Georgia Institute of Tech-nology, 1987.

[124] J. N. Ridley. Computer simulation of contact pressure in capitula.Journal of Theoretical Biology, 95:1–11, 1982.

[125] J. N. Ridley. Packing efficiency in sunflower heads. MathematicalBiosciences, 58:129–139, 1982.

[126] D. F. Robinson. A notation for the growth of inflorescences. NewPhytologist, 103:587–596, 1986.

[127] G. Rozenberg and A. Salomaa. The mathematical theory of L-systems. Academic Press, New York, 1980.

[128] A. Salomaa. Formal languages. Academic Press, New York, 1973.

[129] F. W. Sears, M. W. Zemansky, and H. D. Young. College physics.Addison-Wesley Publ. Co., Reading, 6th edition, 1985.

[130] M. Shebell. Modeling branching plants using attribute L-systems.Master’s thesis, Worcester Polytechnic Institute, 1986.

[131] P. L. J. Siero, G. Rozenberg, and A. Lindenmayer. Cell divisionpatterns: Syntactical description and implementation. ComputerGraphics and Image Processing, 18:329–346, 1982.

222 Bibliography

[132] W. Sierpinski. Sur une courbe dont tout point est un pointde ramification. Comptes Rendus hebdomadaires des seances del’Academie des Sciences, 160:302–305, 1915. Reprinted in W.Sierpinski, Oeuvres choisies, S. Hartman et al., editors, pages99–106, PWN – Editions Scientifiques de Pologne, Warsaw, 1975.

[133] G. Siromoney and R. Siromoney. Rosenfeld’s cycle grammars andkolam. In H. Ehrig, M. Nagl, A. Rosenfeld, and G. Rozenberg,editors, Graph grammars and their application to computer sci-ence; Third International Workshop, Lecture Notes in ComputerScience 291, pages 564–579. Springer-Verlag, Berlin, 1987.

[134] G. Siromoney, R. Siromoney, and T. Robinsin. Kambi kolamand cycle grammars. In R. Narasimhan, editor, A perspective intheoretical computer science — commemorative volume for GiftSiromoney, Series in Computer Science Vol. 16, pages 267–300.World Scientific, Singapore, 1989.

[135] R. Siromoney and K. G. Subramanian. Space-filling curves andinfinite graphs. In H. Ehrig, M. Nagl, and G. Rozenberg, editors,Graph grammars and their application to computer science; Sec-ond International Workshop, Lecture Notes in Computer Science153, pages 380–391. Springer-Verlag, Berlin, 1983.

[136] A. R. Smith. Plants, fractals, and formal languages. Proceedingsof SIGGRAPH ’84 (Minneapolis, Minnesota, July 22-27, 1984)in Computer Graphics, 18, 3 (July 1984), pages 1–10, ACM SIG-GRAPH, New York, 1984.

[137] A. R. Smith. About the cover: Reconfigurable machines. Com-puter, 11(7):3–4, 1978.

[138] P. S. Stevens. Patterns in nature. Little, Brown and Co., Boston,1974.

[139] R. J. Stevens, A. F. Lehar, and F. H. Perston. Manipulationand presentation of multidimensional image data using the Peanoscan. IEEE Trans. on Pattern Analysis and Machine Intelligence,PAMI-5(5):520–526, 1983.

[140] A. L. Szilard. Growth functions of Lindenmayer systems. Techni-cal Report 4, Computer Science Department, University of West-ern Ontario, 1971.

[141] A. L. Szilard and R. E. Quinton. An interpretation for DOLsystems by computer graphics. The Science Terrapin, 4:8–13,1979.

Bibliography 223

[142] R. Thom. Structural stability and morphogenesis. An outline ofa general theory of models. Benjamin/Cummings, Reading, Mas-sachusetts, 1975.

[143] d’Arcy Thompson. On growth and form. University Press, Cam-bridge, 1952.

[144] W. Troll. Die Infloreszenzen, volume I. Gustav Fischer Verlag,Stuttgart, 1964.

[145] W. Troll. Die Infloreszenzen, volume II. Gustav Fischer Verlag,Jena, 1969.

[146] A. Turing. On computable numbers with an application to theEntscheidungsproblem, 1936. Proc. Lond. Math. Soc. (ser. 2),42:230–265, 1936–37, and 43:544–546, 1937.

[147] A. Turing. The chemical basis of morphogenesis. PhilosophicalTrans. Roy. Soc. B, 237(32):5–72, 1952.

[148] W. T. Tutte. Graph theory. Addison-Wesley, Reading, Mas-sachusetts, 1982.

[149] S. Ulam. Patterns of growth of figures: Mathematical aspects. InG. Kepes, editor, Module, Proportion, Symmetry, Rhythm, pages64–74. Braziller, New York, 1966.

[150] J. A. M. van den Biggelaar. Development of dorsoventral polarityand mesentoblast determination in Patella vulgata. Journal ofMorphology, 154:157–186, 1977.

[151] A. H. Veen and A. Lindenmayer. Diffusion mechanism for phyl-lotaxis: Theoretical physico-chemical and computer study. PlantPhysiology, 60:127–139, 1977.

[152] X. G. Viennot, G. Eyrolles, N. Janey, and D. Arques. Combinato-rial analysis of ramified patterns and computer imagery of trees.Proceedings of SIGGRAPH ’89 (Boston, Mass., July 31-August4, 1989), in Computer Graphics 23,4 (August 1989), pages 31–40,ACM SIGGRAPH, New York, 1989.

[153] P. M. B. Vitanyi. Development, growth and time. In G. Rozen-berg and A. Salomaa, editors, The Book of L, pages 431–444.Springer-Verlag, Berlin, 1986.

[154] H. Vogel. A better way to construct the sunflower head. Mathe-matical Biosciences, 44:179–189, 1979.

[155] H. von Koch. Une methode geometrique elementaire pour l’etudede certaines questions de la theorie des courbes planes. Actamathematica, 30:145–174, 1905.

224 Bibliography

[156] J. von Neumann. Theory of self-reproducing automata. Universityof Illinois Press, Urbana, 1966. Edited by A. W. Burks.

[157] F. Weberling. Typology of inflorescences. J. Linn. Soc. (Bot.),59(378):215–222, 1965.

[158] F. Weberling. Morphologie der Bluten und der Blutenstande. Ver-lag Eugen Ulmer, Stuttgart, 1981.

[159] H. Weyl. Symmetry. Princeton University Press, Princeton, NewJersey, 1982.

[160] S. Wolfram. Computer software in science and mathematics. Sci-entific American, 251(3):188–203, 1984.

[161] S. Wolfram. Some recent results and questions about cellularautomata. In J. Demongeot, E. Goles, and M. Tchuente, edi-tors, Dynamical systems and cellular automata, pages 153–167.Academic Press, London, 1985.

[162] T. Yokomori. Stochastic characterizations of EOL languages. In-formation and Control, 45:26–33, 1980.

[163] D. A. Young. On the diffusion theory of phyllotaxis. Journal ofTheoretical Biology, 71:421–423, 1978.

[164] M. H. Zimmerman and C. L. Brown. Trees — structure andfunction. Springer-Verlag, Berlin, 1971.

Index

accumulation 67acrotonic structure 87age 135alphabet 4Anabaena catenula 5, 136anastomosis 51angle increment 6apex 21

continuing 70flowering 64terminal 70vegetative 64

apical dominance 88apical front 155apple 75attractor 178axiom 4axis 21

bicubic patch 119blastula 166blue-green bacteria 5branch 21branching

decussate 93monopodial 55, 70polypodial 70, 86sympodial 50, 58, 70ternary 58

branching angles 52branch plane 52

capitulum 96daisy 102sunflower 102

Capsella bursa-pastoris 73, 120Carex laevigata 116catastrophe 134cell 150

cell (continued)apical 154epidermal 173marginal 163vegetative 44

cell division system 154cellular automata 51cellwork 168chaos game 178cleavage 166condition 42contact point 13, 119context 30continuity requirement 140control graph 182crocus 68curve

dragon 11E-curve 12FASS 12Gosper 12Hilbert 16, 21Koch island 7Peano 17pure 16Sierpinski 11snowflake 1space-filling 12SquaRecurve 14

cyme 82closed 85double 83open 82simple 82

data base amplification 63da Vinci’s postulate 57delay mechanism 66derivation 4, 42, 135, 147, 171

225

226 Index

derivation (continued)direct 4length 4stochastic 28

developmental sequence 4developmental switch 65dibotryoid 76

closed 79open 76

direction vector 13divergence angle 52, 100, 109Dryopteris thelypteris 162

edge 168directed 146neutral 146

embryo 166encyclic number 110entry point 13entry vector 13Euler method 153exit point 13exit vector 13expression 41

fern 120Fibonacci

angle 99number 99, 102series 37

florigen 76flower head 96flowering sequence

acropetal 73, 78basipetal 75, 78, 87

fractal 6

gametophyte 153generative helix 109generator 1golden mean 99grammar 2growth

exponential 36, 142linear 140polynomial 37, 143sigmoidal 38, 39square-root 38

growth function 36, 140growth potential 88, 123

heading 6heading vector 119heterocyst 44Hooke’s law 151, 172

IFS 177controlled 182

inference 11inflorescence 71initiator 1interaction 65internode 21

Koch construction 1, 7

L-system 21L 302L 30bracketed 24cellwork 168complete 64context-sensitive 30deterministic 4DOL 4IL 30map 146mBPCOL 168mBPMOL 146, 168OL 4parametric 41partial 64pseudo 18scheme 64stochastic 28, 66string 4table 66timed 135tree 23

leafcompound 128, 142cordate 122simple 123

letter 3timed 135

letter occurrence 28

Index 227

lifetime 135lilac 93lily-of-the-valley 73, 127lineage 65Lychnis coronaria 83

map 145map interpretation

center of gravity 150dynamic 151, 172wall subdivision 150

marker 147, 170matching 147

Microsorium linguaeforme 153microworld 187mint 79module 41, 134Mycelis muralis 87

Newton’s law 152nitrogen compound 44

object 189offset 148order of continuity 141osmotic pressure 151, 173

panicle 86parameter 41parastichy 102

order 110Pascal triangle 38Patella vulgata 166periclinal ratio 166phase effect 63, 73phyllotaxis 99

cylindrical 109planar 100

pineapple 116plastochron 65predecessor 4, 42, 147production 1

edge 146, 168identity 4inconsistent 171probabilistic 28string 4tree 23

raceme 71closed 75compound 76, 81open 71

recursive formula 15region 145rewriting 1

edge 11, 46node 13, 47

root node 21rose 106, 126

segment 154basal 159

signal 31, 67acropetal 76, 88basipetal 88

spadix 96spike 95spruce cone 116statement

#define 44#ignore 32#include 44

step size 6string 3

bracketed 24subapical growth 26subfigure 13successor 4, 42, 147

terminal node 21thallus 154thyrsus 85tile 12total elapsed time 136tree

axial 21rooted 21

tribotryoid 81tropism vector 58turtle interpretation 6, 18, 26,

46, 57, 73, 119, 120, 122,127

umbel 95Umbelliferae 128up vector 119

228 Index

virtual laboratory 187voxel space 51

wall 150, 168anticlinal 155periclinal 155tension 151, 173

water-lily 106wild carrot 95, 131word 3

parametric 41timed 135

zinnia 106