# Riemannian L-systems:
Modeling growing forms in curved spaces

Christophe Godin  Laboratoire Reproduction et Développement des Plantes, Univ
Lyon, ENS de Lyon, UCB Lyon1, CNRS, INRAE, Inria, F-69342 Lyon, France
Frédéric Boudon  CIRAD, UMR AGAP Institute, F-34398 Montpellier, France  UMR
AGAP Institute, Univ. Montpellier, CIRAD, INRAE, Institute Agro, Montpellier,
France

###### Abstract

In the past 50 years, the formalism of L-systems has been successfully used
and developed to model the growth of filamentous and branching biological
forms. These simulations take place in classical 2-D or 3-D Euclidean spaces.
However, various biological forms actually grow in curved, non-Euclidean,
spaces. This is for example the case of vein networks growing within curved
leaf blades, of unicellular filaments, such as pollen tubes, growing on curved
surfaces to fertilize distant ovules, of teeth patterns growing on folded
epithelia of animals, of diffusion of chemical or mechanical signals at the
surface of plant or animal tissues, etc. To model these forms growing in
curved spaces, we thus extended the formalism of L-systems to non-Euclidean
spaces. In a first step we show that this extension can be carried out by
integrating concepts of differential geometry in the notion of turtle
geometry. We then illustrate how this extension can be applied to model and
program the development of both mathematical and biological forms on curved
surfaces embedded in our Euclidean space. We provide various examples applied
to plant development. We finally show that this approach can be extended to
more abstract spaces, called abstract Riemannian spaces, that are not embedded
into any higher space, while being intrinsically curved. We suggest that this
abstract extension can be used to provide a new approach for effective
modeling of tropism phenomena and illustrate this idea on a few conceptual
examples.

##  1 Introduction

Morphogenesis is the process by which biological and non-biological forms
develop over time. Systems whose state changes in time according to some
evolution rules are called dynamical systems, eg. see [72]. These rules are
usually formalized as a differential equation, called the evolution equation,
that expresses the rate of change of the system’s state as a function of the
current system’s state and time:

$\frac{\partial\mathbf{x}(t)}{\partial
t}=F_{\lambda}(\mathbf{x}(t),\partial_{i}\mathbf{x}(t),t,\mathbf{z}(t)),$ (1)

where $\mathbf{x}(t)$ denotes the system’s state, usually in $\mathbb{R}^{n}$
for some integer $n$, $\partial_{i}\mathbf{x}(t)$ denotes its spatial
variation rate in the $i=1\cdots n$ coordinate directions of space,
$F_{\lambda}$ specifies the rule describing how the system in state
$\mathbf{x}(t)$ changes over a small time step $dt$, $\lambda$ refers to the
model parameters, and $\mathbf{z}(t)$ denotes any external variable that may
affect the dynamics of the system.

When the system is a form that evolves over time, the state $\mathbf{x}$
represents the form as a mathematical structure (e.g. a sequence, a tree, a
graph, a mesh, …) loaded with additional state variables at the different
positions in the structure, and forming so called fields on the form [19].
Such dynamic systems have a state whose structure changes over time and are
called dynamical systems with dynamical structures [17]. In this case, Eq.1
must be generalized to account for the more complex nature of the system’s
state and the lack of a notion of derivative on discrete structures:

$\mathbf{x}(t+dt)=F_{\lambda}(\mathbf{x}(t),\partial_{i}\mathbf{x}(t),t,\mathbf{z}(t),dt),$
(2)

where $\mathbf{x}(t)$ denotes this time the system’s structure and its fields,
and $dt$ a small time increment.

Modeling the evolution of a form over time thus comes down to finding a
function $F_{\lambda}$ able to describe the variations of the structure
$\mathbf{x}(t)$ over time. Despite promising efforts [18, 44], there is yet no
general theory to define expressions of $F_{\lambda}$ for general form growth.
However, restrictions to specific form structures or dimensions have lead to
the development of efficient formalisms to model form development, see [61,
21] for reviews. This is particularly well illustrated by the formalism of
L-systems, [37, 38, 39, 59], that was introduced to model the development of
sequences and branching structures, very common in biology.

L-systems have been used successfully over more than 50 years to model the
construction of fractal forms and the growth of plant branching structures,
such as plant architectures, inflorescences, vein patterns, root systems, etc.
L-systems make it possible to model the evolution of both the structure and
the geometry of forms. For this, the state of a growing form, i.e. a sequence
or branching structure augmented with field variables, is mapped to a 3-D
shape using turtle geometry [1]. Remarkably, the formalism leads to the
definition of a computer language that can be naturally used to program the
development of forms in a declarative manner. Various computational
implementation of this formalism have been proposed, e.g. [59, 57, 25, 6], as
extensions of different programming languages, and putting emphasis on the
development of different aspects of the theory.

Common to all these implementations is the assumption that the geometric
interpretation of the forms is carried out in an Euclidean (flat) space as a
vast majority of plant development models are conceived in flat spaces. Most
models of plant architecture development for instance simulate the growth of a
branching structure in 3-dimensional (3-D) Euclidean space, e.g. [60, 20, 51,
5]. Models of vein development in leaves simulate the growth of a branching
vascular structure within 2D shapes representing a flat leaf blade, e.g. [13,
42, 67]. In plant tissues, models of hormone signaling or water flows usually
assume transport laws expressed in 2-D or 3-D Euclidean representations of the
tissue, e.g. [30, 24, 71, 8]. Likewise, the simulation of filamentous system
growth, such as pollen tubes or root hairs, is carried out essentially within
flat embedding substrate spaces, e.g. [12, 10].

Yet, a number of these phenomena actually take place in non-flat and so-called
curved spaces. Climbing plants for instance may growth on tree trunks or
grounds that are not flat surfaces. Vein networks can develop in leaf blades
that are markedly curved. Pollen tubes grow on the pin-like structures of
papillae that are not flat [64]. Microtubules polymerize/depolymerize
dynamically within the cell cortex that in general is not flat [2].

Examples are numerous and occur on a variety of scales. However, the modeling
of plant form growth in curved spaces has been up to now only scarcely
investigated. At the level of macro-molecules, simulations of microtubule
dynamics have been carried out in 3D cell geometries to study the emerging
properties of such a network of microfilaments subjected to local synthesis,
decay and interaction rules, [43]. In this work, cell geometry is represented
by a 3-D mesh, and microtubule trajectories are computed by assuming that
microtubules are progressing in straight line in the 3-D Euclidean space. The
resulting displacement is then projected the local tangent plane to account
for potential curvature of the cell geometry. Likewise, at organ scale a
similar projection strategy is used in [26] to model climbing plants or in
[65] to model branching venation patterns on the surface of petals. In both
cases, to approximate the formation of a pattern in a curved space, the
pattern growth is first evaluated in the 3-D Euclidean ambient space. The
resulting form is then projected on the discrete curved surface represented as
a 3-D triangular mesh. A more direct use of the concepts of differential
geometry was described in a different application context to create artistic
patterns on the surface of 3D objects by [36]. This approach exploits user-
defined vector fields on surface meshes to locally drive drawings of curves at
the surface. This is different from the approach that we use here which is
based on the possibility to follow geodesics in curved spaces to construct
forms. However, similarly to what we do here, the authors construct a language
that makes it possible to program patterns based on vector fields. A different
approach, aimed at modeling the growth of plant lianas of their support and
more generally the growth of parasite organisms on various types of hosts,
analyzes how to grow surfaces by accretion formalized in [46] with the
constraint of keeping on a reference surface representing the host [49]. In
this approach, the trajectories of the parasite is defined analytically on the
host surface and the focus is on the construction of the surface representing
the parasite envelop. By contrast, our approach is primarily focused on the
construction of the trajectories representing the patterns growing in curved
spaces. In a preliminary work, we explored the possibility of using L-systems
to model fractal structures on simple spheres [62]. Here, we largely extend
this initial exploration with a complete formalization of the notion of
Riemannian L-systems, that can be applied to both smooth curved surfaces and
more abstract non-Euclidean smooth curved spaces.

Our approach couples L-systems and differential geometry. Importantly, this
extension remains easy to use for modelers as it allows to program L-system
models as if processes were locally taking place in a Euclidean space. To draw
geometric patterns, the user mainly thinks in terms of elementary movements,
’go straight’, ’turn right’, etc. without paying attention most of the time to
the curvature of the underlying space. In this way, we show that turtle
geometry makes it possible to define a notion of intrinsic shape, that does
not depend on the embedding space. Specific primitives allow the user to use
geometric properties of curved spaces (curvature, excess angle, parallel
transport, point-wise geodesic distance, etc.) to develop programs in which
the geometry of the embedding space feeds back on the developing form. Through
the paper, we chose not to assume that the reader is familiar with concepts in
differential geometry. We therefore introduce the necessary fundamental
concepts and notations used in this domain to keep the text as self-contained
as possible. We also provide various applications of how Riemannian L-systems
can be used to illustrate key concepts of differential geometry and to explore
of a variety of mathematical or biological dynamical systems in curved spaces
such as the growing fractal forms, random walks, developing branching
structures, tip growing filaments, vein pattern generation and so on.

##  2 L-systems overview

L-systems are basically rewriting systems on strings for which rewriting rules
are applied in parallel to all the elements of the current string to compute
the new string. In this section, we briefly introduce L-systems [39, 59] and
key associated notions.

###  2.1 Basic formalism

Let $V=\\{a_{1},a_{2},a_{3},...,a_{N}\\}$ be a finite set of elements (called
symbols or modules). We call $V^{*}$ the set of all finite strings that can be
constructed by concatenating any number of symbols from $V$ ($V^{*}$ is the
free monoid constructed from $V$ for the binary operation of string
concatenation). $V^{*}$ includes the empty string, denoted $\lambda$, that is
the neutral element for string concatenation. Sequences of $V^{*}$ are called
words. For example, if $V=\\{A,B,a,b\\}$, $x_{1}=aaBAa$, and
$x_{2}=abbAaAabbbb$ are words of $V^{*}$.

######  Definition 1 (String homomorphism).

Let us consider a word $x\in V^{*}$ and a decomposition of this word into
subwords: $x=x_{1}x_{2}\cdots x_{K}$, where $x_{k}$ are also words $\in
V^{*}$, for $k=1\dots K$. A homomorphism $H$ from $V^{*}$ to $V^{*}$ is a
mapping such that:

$H(x_{1}x_{2}\dots x_{K})=H(x_{1})H(x_{2})\dots H(x_{K}).$

From this definition, it follows that a string homomorphism is completely
defined by the definition of the images of the symbols in $V$. In addition,
this definition implies that $H(\lambda)=\lambda$.

Let us call $w_{n}$ the word image by $H$ of symbol $a_{n}\in V$,
$H(a_{n})=w_{n}$. In the context of L-systems, the pair $P_{n}=(a_{n},w_{n})$
is called a production rule. We denote $P=\\{P_{n}\\}_{n=1\cdots N}$ the set
of all production rules associated with $H$.

######  Definition 2 (D0L-system).

A D0L-system $\mathcal{L}$ is a 3-uple $(V,P,A)$, where
$V=\\{a_{1},a_{2},a_{3},...,a_{N}\\}$ is a finite set of symbols, $P$ is a set
of production rules on $V$, and $A\in V^{*}$ is called the axiom.

######  Definition 3 (Derivation).

Let $V=\\{a_{1},a_{2},a_{3},\cdots,a_{N}\\}$ be a finite set of symbols, $P$ a
set of production rules on $V$, and $H$ the homomorphism associated with $P$.
Let $x$ be a word in $V^{*}$, $H(x)$ is called the derivation of $x$ by $P$.
If $x=x_{1}x_{2}\cdots x_{K}$ and $H(x_{k})=w_{k}$ for $k=1\dots K$, the
derivation of $x$ is denoted by:

$x_{1}x_{2}\cdots x_{K}\Rightarrow w_{1}w_{2}\cdots w_{K}.$

For a D0L-system $\mathcal{L}=(V,P,A_{0})$, successive derivations of the
axiom $A_{0}$ can be obtained in a deterministic manner. Let us denote $A_{i}$
the $i$-th derivation of $A$, i.e. $A_{i}=H^{i}(A_{0})$, we have:

$A_{0}\Rightarrow A_{1}\Rightarrow A_{2}\Rightarrow\cdots\Rightarrow A_{i}.$

A (D0)L-system may thus represent dynamical systems whose states can be
abstracted as (discrete) strings of components, called L-strings, that
dynamically change as derivations are applied. Derivations are often
interpreted as the evolution of the system’s state with time. They can also
represent changes of the observer’s viewpoint (such as zooming in a structure,
which shows more elements, and thus result in a change of the system’s
representation).

##### Example

A simple example of a D0Lsystem is provided by the development of a
multicellular filamentous organism Anabenae [59] based on a model originally
developed by Koster and Lindenmayer [33]. This organism is organized as a file
of cells that have different types ${a,b,A,B}$ corresponding to their
differentiation and polarization states. According to their types, cells have
different developmental behaviors described by specific rules that can be
modeled by a L-system [59].

Let $\mathcal{L}=(V,P,A_{0})$, where $V=\\{A,B,a,b\\}$, $P=\\{A\rightarrow
Ba,B\rightarrow bA,a\rightarrow A,b\rightarrow B\\}$, and $A_{0}=A$. Starting
from the axiom $A_{0}=A$, the sequence of derivations of this L-system goes as
follows:

$A\Rightarrow Ba\Rightarrow bAA\Rightarrow BBaBa\Rightarrow
bAbAAbAA\Rightarrow\cdots$

Each string (L-string) of this derivation sequence represents a state of the
growing organism at consecutive time steps.

###  2.2 D0L-system extensions

D0L-systems have been extended in various ways and have lead to develop
powerful languages constructs to simulate dynamical systems. Here we briefly
recapitulate common key extensions [59].

##### Branching systems.

First, D0L-systems have been extended to model branching systems rewriting and
not only strings. This extension relies on the fact that branching systems can
be simply encoded as bracketed strings. This makes it possible to define
readily D0L-systems that rewrite branching systems. For this, the definition
of the vocabulary $V$ is augmented by a pair of opening ’[’ and closing ’]’
square brackets. In addition, a restriction is imposed on the homomorphism on
$V^{*}$: square brackets can only be mapped to themselves by $H$, and they can
only appear on the right-hand side of a production rule if they form a well-
formed bracketed string (all opening brackets must be properly balanced in the
string and opening/closing brackets must be strictly nested).

In a plant for instance the apex A of a stem can produce a new portion of stem
(internode I), a lateral apex A and a new apical apex. This can be represented
by the production rule:

$A\rightarrow I[A]A.$

Starting from the axiom $A$, the first derivations yield:

$A\Rightarrow I[A]A\Rightarrow I[I[A]A]I[A]A\Rightarrow\cdots$

Note that it is usually assumed that every module for which a production rule
is not specified (here $I$ for example) is rewritten unchanged in the new
string (identity transformation).

##### Parametric rules.

A second powerful extension is the possibility to add parameters to the
modules. This makes it possible to write rules that propagate parameter values
as the modules are rewritten and to make computation on them. For instance,
the elongation of a rectangular cell represented by a module $C$ can be
modeled by introducing a real parameter $y$ in a production rule such that:

$C(y)\rightarrow C(y+\delta y),$

where $\delta y$ is an increment of length. If the axiom consists of the
string $C(0)$, then applying the above production rule to the axiom will yield
the following $i$ first derivations:

$C(0)\Rightarrow C(\delta y)\Rightarrow C(2\delta
y)\Rightarrow\cdots\Rightarrow C(i\delta y).$

##### Context-sensitive rules.

Another useful extension is the notion of context-sensitive rules. Here, a
production rule is applied to the left-hand side module only if the context of
this module matches some criterion in the current L-string. Traditionally, the
left and right context of a module in a string are specified by using ’<’ and
’>’ markers.

For instance, the rule:

$B<A\rightarrow B,$

means that a symbol A must be rewritten into B only if its left-context (i.e.
the module immediately to its left) in the original L-string is a B.
Otherwise, the symbol is left unchanged. The action of this rule can be
observed on the axiom $BAA[AA]AA$ that leads to the sequence of derivations:

$\displaystyle BAA[AA]AAA$ $\displaystyle\Rightarrow BBA[AA]AAA\Rightarrow
BBB[AA]AAA$ $\displaystyle\Rightarrow BBB[BA]BAA\Rightarrow
BBB[BB]BBA\Rightarrow BBB[BB]BBB.$

This makes it possible for instance to model the propagation of a signal in a
branching structure. Such a signal can propagate from the root of the
branching structure to the leaves using left-context rules, or from the leaves
to the root using right-context rules. One can note that the rule
$B<A\Rightarrow B$ applies to L-string patterns such as $..BA..$ or $..B[A]..$
and thus take into account the branching organization of the L-string.

##### Procedural rules.

Instead of specifying once for all the right-hand side of production rules, it
might be interesting to determined it procedurally at runtime. For example,
the procedural production rule:

$\displaystyle A(x):$ $\displaystyle\mathrm{if}\;x>0.3:$
$\displaystyle\rightarrow B,$ $\displaystyle\mathrm{else}:$
$\displaystyle\rightarrow C.$

This rule yields different derivations for different axioms:

$\displaystyle A(0.5)$ $\displaystyle\Rightarrow B,$ $\displaystyle A(0.1)$
$\displaystyle\Rightarrow C.$

Procedural rules thus have a classical left hand-side, and a right-hand side
that is a procedure ending with rewriting statements (here $\rightarrow B$ and
$\rightarrow C$). Procedural rules make it possible for instance to simulate
non-deterministic L-systems. In this case the right-hand side rewriting
statement is computed based on some random choice.

All the above extensions can be combined to model complex dynamical
behaviours. Note that the resulting mathematical structures are no longer
D0L-systems, but extensions of them that will be called hereafter using the
generic term of L-systems.

Based on the pioneering language called cpfg developed by P. Prusinkiewicz
[59] a number of language variants have been built over the years on the top
of different programming languages, e.g. L+C (C++) [57], XL (Java)[32], L-Py
(Python)[6]. In this paper, the examples and extensions are developed in L-Py.

###  2.3 Turtle geometry: adding geometry in L-systems

L-strings are abstract representations of a system’s state with no particular
graphic representation. It is however often very useful to attach a geometric
representation of the system’s state in the 3D-space (to display the 3D
architecture of a modeled growing plant for example). For this L-systems make
use of turtle geometry [1].

##### Turtle geometry

Turtle geometry is a way to define complex geometric objects in 3D as a
sequence of elementary geometric instructions. Basically, a (virtual) turtle
is able to move and draw in the 3D space as it moves. For this, the turtle is
given a sequence of elementary geometric instructions that it can read and
interpret sequentially. Interestingly, as the definition of forms using turtle
geometry relies on purely geometric primitives, the construction process is in
general independent the selected coordinate system.

We assume that the turtle moves with respect to a global reference frame,
denoted as $\mathcal{R}_{0}$ and that it is itself represented as a moving
local reference frame, hereafter called the turtle frame
$\mathcal{R}=\\{\mathbf{H},\mathbf{L},\mathbf{U}\\}$ (respectively denoting
unit vectors related to the turtle’s body: Head, Left, Up). Each elementary
instruction is interpreted by the turtle as an order to move or to draw [1,
59]. The elementary instructions are coded as strings of geometric modules
having specific names. For example, the module F(l) instructs the turtle to
move forward by a distance $l$ in the direction $\mathbf{H}$ of its head, and
draw a line while moving, the module +(a) means that the turtle should turn
left (around its Up axis $\mathbf{U}$) by a angle $a$ degrees, the module ;(c)
means that color $c$ should be used for drawing from now on, etc.. A L-system
language usually contains a number of predefined geometric modules that make
it possible to draw a large variety of simple and more complex geometrical
shapes (see detailed list of turtle instructions in L-Py for example).

To operate on an input string, the turtle is associated with a state $S$ that
records its current information: position $(x,y,z)$, orientation
$\\{\mathbf{H},\mathbf{L},\mathbf{U}\\}$ expressed in the global reference
frame $\mathcal{R}_{0}$, its current color $c$ , etc.,

$S=(x,y,z,H,L,U,c,\cdots).$

When reading a new instruction from an input sequence of turtle commands, the
turtle executes the elementary action corresponding to the read string module,
updates its state accordingly and proceeds to the next instruction in the
sequence. Sequences of geometric modules, interpretable by a turtle, are
called T-strings.

##### Branching systems.

As explained above, branching systems can be modeled by using square brackets.
During turtle interpretations, when reaching an opening square bracket in an
L-string, the turtle saves its current state on the top of a stack, called the
interpretation stack, and then proceeds with the string interpretation inside
the bracket. When reaching a closing square bracket, the turtle pops the
current top state of the interpretation stack and restores its current state
with it before continuing to interpret the L-string. This push/pop turtle
mechanism linked with the use of square brackets makes it possible to easily
create branching patterns in L-systems. To insert a branch on another for
example, a well balanced bracketed list of modules representing the new branch
must be inserted at the position corresponding to the bifurcation between a
main branch and the new lateral branch.

##### Interpretation rules: Coupling L-strings with turtle geometry

L-systems strings can be associated with a geometry using turtle geometry.
This can be done by either directly integrating symbols that can be
interpreted by the turtle in the finite set of symbols $V$ or by using
interpretation rules to translate L-strings (corresponding to the dynamic
system’s state) into T-strings (corresponding to its geometrical
interpretation).

Remarkably, these geometric interpretation rules can also be represented as an
homomorphism $G$ mapping L-strings to T-strings [34]. This homomorphism is
usually defined on set of L-system symbols $V$. Let us call $z_{n}$ the word
image by $G$ of symbol $b_{n}\in V$, $G(b_{n})=z_{n}$. The pair
$I_{n}=(b_{n},z_{n})$ is called a interpretation rule. We denote
$I=\\{I_{n}\\}_{n=1\cdots M}$ the set of all the interpretation rules
associated with $G$.

The interpretation rules are thus used to make a translation between the
L-string modules that in general denote the components of the modeled system
with their attributes and have no direct geometric interpretation, and the
T-string geometric modules that can be directly interpreted by the turtle as
geometrical instructions.

Interpretation rules can be recursive and are applied in a depth-first manner
to the input L-string (recursions are fully developed before processing to the
next module of the string). Most of the extensions that have been defined for
L-system derivation rule above are also defined for interpretation rules
(parametric, procedural), except for context sensitivity.

##### Sensing the environment.

In L-systems, forms get developed in a 2D or 3D Euclidean embedding space. As
they grow, it may sometimes be useful to locally probe the environment in a
production rule to make a decision that will impact the growth. For example,
one may want to determine if a new tentative segment extension of the growing
form collides with some object in the environment or with previously
constructed part of the same form to decide if the grow will actually take
place.

For this, it is essential to access the position and orientation of the turtle
in production rules. A specific mechanism has thus been designed in this aim
[56]. It relies on the use of specific query modules such as ?P ?H or ?U, that
allow the modeler to access the turtle position and orientation at different
locations of a L-string from within production rules. To access such
information at a given derivation step, these query modules have to be
produced in the L-string at the previous iteration step. They act as place-
holders to record geometric information computed during turtle’s
interpretation at the previous step. Once filled, the parameters of these
modules can be used in production rules of the next iteration to make
development decisions.

###  2.4 L-system examples in L-Py

Let us illustrate the above concepts, on a few examples showing how simple
forms can be computed in space or in both space and time with L-systems,
Fig.1, together with their L-system code (here using L-Py).

The first example illustrates the notion of recursive production rule, where a
module A containing a variable $n$ produces a segment of size $n$, turns of a
fixed angle and produces a new module A with variable $n+1$. When repeated
(here 20 times) this produces a polygonal Archimedean spiral, Fig.1a. On the
same principle, the second example shows how stochasticity can be introduced
in the rules to generate a random walk for instance 1.b. The use of
recursivity is at the core of the language. If extended beyond spatial
recursivity, it can be used to model a wide variety of shapes. In Figure1.c,
for example, recursivity is used over scales to produce a fractal form [41].
Here, at each derivation, existing segments are recursively rewritten as
several segments of shorter size. In this procedure the turtle interpretation
of the Lstring converges towards a fractal form [59]. Finally, the recursive
principle can also be used in space and time to simulate the development of a
form. Figure1.d illustrates how the development of a plant branching system
can be described very concisely using this principle.

⬇

1Axiom: A(1)

2derivation length: 20

3production:

4A(n): nproduce F(n)+(60)A(n+1)

Listing 1: Archimedean spiral (see Fig.1a)

⬇

1dl = 1.

2Axiom: A(0)

3derivation length: 1000

4production:

5A(n):

6 a = 360*random()

7 nproduce +(a)F(dl)A(n+1)

Listing 2: Random walk in 2-D (see Fig.1b)

⬇

1Axiom: F(1)-(120)F(1)-(120)F(1)

2derivation length: 5

3production:

4F(x) : nproduce F(x/3.0)+F(x/3.0)--F(x/3.0)+F(x/3.0)

Listing 3: Fractal curve (von Koch flake) (see Fig.1c)

⬇

1N = 10

2offset = 2

3iangle = 60

4

5Axiom: A(0)

6derivation length: N

7production:

8A(n) :

9 if n<N: # produces an internode, a lateral bud (in [...]) and an apical bud

10 nproduce I(n) [P(n)A(n+offset)] A(n+1)

11 else: # produces a flower bud

12 nproduce B

13

14interpretation:

15I(n) : nproduce ;(1)F(N-n)

16P(n) : # Phyllotaxis angle

17 if n%2 == 0 : nproduce +(iangle)

18 else: nproduce -(iangle)

19A(n): nproduce ;(2)@O(1.5)

20B : nproduce F(1);(3)@O(1.5)

Listing 4: Plant branching structure development (see Fig.1d) Figure 1:
Examples of geometric forms that can be produced by Lsystems (here using the
computer language LPy). (a) a polygonal Archimedean spiral, (b) a random walk,
(c) the von Koch flake (fractal curve) (d) an idealized simple sequence of
development of a plant branching system.

Forms built using classical L-systems assume that the underlying 2D or 3D
space is Euclidean, i.e. that the turtle graphical commands are interpreted as
if the turtle were moving in a Euclidean space. Euclidean spaces are flat,
i.e. the five Euclid postulates that found classical geometry hold. In
particular, the fifth one, called the parallel postulate, that states that, in
a plane, through a point exterior to a given straight line, at most one line
passes that never intersects the initial line. It can be shown that flat
spaces are characterized by this fifth postulate [48]. In the sequel, we
extend L-systems to operate in curved (non flat) spaces and provide a
programming language to describe form and form development in these more
general non-Euclidean spaces.

##  3 Moving on parametric surfaces

In this aim, we start by studying the extension of L-systems to 2-D curved
surfaces embedded in a 3-D Euclidean space. For this, we explore how turtle
geometry, which defines how forms are constructed in space, can be extended to
2-D surfaces. The extension of turtle geometry to non-parametric, mesh-like
surfaces has been described in [1]. Here we investigate the extension of
turtle geometry to parametric surfaces.

###  3.1 Parametric surfaces

A (smooth) parametric surface $\mathcal{S}$ is defined as a differentiable map
from a 2-D parameter space $\mathbb{U}\subset\mathbb{R}^{2}$ to
$\mathbb{R}^{3}$. Let us denote $(u^{1},u^{2})$ the coordinates of points in
$\mathbb{U}$ and $x^{1},x^{2},x^{3}$ the coordinates of points in
$\mathbb{R}^{3}$. With these notations, a 2-D surface can be defined by the
equations:

$\begin{array}[]{ll}x^{1}=&\phi^{1}(u^{1},u^{2})\\\
x^{2}=&\phi^{2}(u^{1},u^{2})\\\ x^{3}=&\phi^{3}(u^{1},u^{2}).\\\ \end{array}$

These 3 functions are often summarized by writing more simply
$\mathbf{x}(\mathbf{u})$, or in coordinates
$x^{i}(u^{\alpha})=\phi^{i}(u^{\alpha})$, $i=1,2,3$, $\alpha=1,2$ reminding us
that $x^{i}$’s are functions of the $u^{\alpha}$’s. If $u^{1}$ (resp. $u^{2}$)
has a fixed value, the variation of the other parameter $u^{2}$ (resp.
$u^{1}$) defines a so-called coordinate line on the surface, Fig.2a. Any pair
of parameters $(u^{1},u^{2})$ thus defines a point $P$ on the surface. The
spatial infinitesimal variations of this point $P(u^{1},u^{2})$ on the surface
with respect to the parameter coordinates, defines two vectors:

$\mathbf{e}_{\alpha}=\frac{\partial P}{\partial u^{\alpha}},$ (3)

that form the covariant basis at point $P$ (also called the coordinate basis).
These vectors are tangent at point $P$ to the coordinate lines, and form a
basis of the plane $T_{P}\mathcal{S}$ tangent to the surface at point $P$. A
vector $\mathbf{X}$ in $T_{P}\mathcal{S}$ is thus a linear combination of the
coordinate basis vectors $\mathbf{e}_{\alpha}$ at $P$ (Fig.2b), and we note,
using Einstein’s implicit summation convention:

$\mathbf{X}=\sum_{\alpha}X^{\alpha}\mathbf{e}_{\alpha}=X^{\alpha}\mathbf{e}_{\alpha}.$
(4)

A vector field over the surface associates a vector in $T_{P}\mathcal{S}$ with
each point $P$ of the surface. Each vector can thus be decomposed in the local
covariant basis following Equ.4.

Figure 2: Manifold curvilinear coordinates illustrated on a manifold of
dimension $2$ embedded in a Euclidean pace of dimension $3$. (a) Coordinate
lines together with the covariant basis at a point $P$. (b) A vector
$\mathbf{X}$ in the tangent plane $T_{P}\mathcal{S}$ at point $P$ can be
decomposed in the covariant basis:
$\mathbf{X}=X^{\alpha}\mathbf{e}_{\alpha}=X^{1}\mathbf{e}_{1}+X^{2}\mathbf{e}_{2}$.

The embedding Euclidean space, $\mathbb{R}^{3}$, induces a metric on the
surface by assigning a value to the dot product of each pair of vectors in the
tangent plane $T_{P}\mathcal{S}$ at $P$:

$<\mathbf{X},\mathbf{Y}>=<X^{\alpha}\mathbf{e}_{\alpha},Y^{\beta}\mathbf{e}_{\beta}>=X^{\alpha}Y^{\beta}<\mathbf{e}_{\alpha},\mathbf{e}_{\beta}>=X^{\alpha}Y^{\beta}g_{\alpha\beta},$
(5)

where, $<.,.>$ denotes the scalar product in $\mathbb{U}^{2}$ and
$g_{\alpha\beta}$ defines the surface metric tensor as the dot product of
$\mathbb{R}^{3}$ vectors:

$g_{\alpha\beta}=<\mathbf{e}_{\alpha},\mathbf{e}_{\beta}>=\mathbf{e}_{\alpha}.\mathbf{e}_{\beta}.$
(6)

The inverse metric tensor $g^{\alpha\beta}$ is such that
$g_{\alpha\gamma}g^{\gamma\beta}=\delta^{\beta}_{\alpha}$, where
$\delta^{\beta}_{\alpha}$ is the Kronecker symbol (= 1 if $\alpha=\beta$ and 0
otherwise) . At each point $P$ of the surface we can also define a normal unit
vector $\mathbf{n}$, perpendicular to all vectors in the tangent plane
$T_{P}$:

$\mathbf{n}=\frac{\mathbf{e}_{1}\times\mathbf{e}_{2}}{|\mathbf{e}_{1}\times\mathbf{e}_{2}|}.$
(7)

We assume in this paper that the surface is orientable, meaning that the
normal vector at each point can be defined in a unique and consistent manner
throughout the surface (this discards non-orientable surfaces such as the
Moebius strip from our analysis.). The orientation is made so that the three
vectors $(\mathbf{e}_{1},\mathbf{e}_{2},\mathbf{n})$ have a direct
orientation, i.e. $\det(\mathbf{e}_{1},\mathbf{e}_{2},\mathbf{n})>0)$.

Surfaces are in general curved. Let us briefly recall what does this mean.
Consider a point $P$ on a surface $\mathcal{S}$ and the normal vector
$\mathbf{n}$ to the surface at $P$, Fig.3ab. Let us consider a direction
$\mathbf{t}_{\theta}$ at $P$ making an angle $\theta$ with a fixed arbitrary
direction $\mathbf{t}_{0}$ in the tangent plane at $P$, and construct the
’vertical’ plane passing by $\mathbf{t}_{\theta}$ and $\mathbf{n}$ at $P$ (in
grey on Fig.3ab). The intersection of this ’vertical’ plane and the surface is
a curve $\gamma_{\theta}$. We assume that a point moves at constant speed on
this curve parameterized by $s$. The velocity of this point is colinear with
$\mathbf{t}_{\theta}$ at $P$ and its rate of variation, defining the curve
curvature at $P$, is thus a vector perpendicular to $\mathbf{t}_{\theta}$:

$\frac{d\mathbf{t}}{ds}=\mathbf{k}_{\theta}.$ (8)

This curvature vector lies in the ’vertical’ plane which contains the curve
$\gamma_{\theta}$. When the direction $\mathbf{t}_{\theta}$ is rotated around
the vertical axis, the norm $k_{\theta}$ of $\mathbf{k}_{\theta}$ varies
continuously. It can be shown, that it passes by a minimum and a maximum
value, called principal curvatures and denoted $k_{min}$ and $k_{max}$, in
specific directions, $\mathbf{d}_{min}$ and $\mathbf{d}_{max}$ called the
principal curvature directions at $P$. As the surface is oriented, the
principal curvatures $k_{min}$ and $k_{max}$ may be either positive (curvature
vector oriented like $\mathbf{n}$) or negative (curvature vector oriented in
the other direction). Principal curvature directions have the property to
always be perpendicular [22].

From the principal curvatures at $P$, one can define the mean curvature
$\kappa_{M}$ and the Gaussian curvature $\kappa_{G}$ as follows:

$\begin{split}\kappa_{M}=&\frac{1}{2}(k_{min}+k_{max}),\\\
\kappa_{G}=&k_{min}\;k_{max}.\end{split}$ (9)

The curved geometry of surfaces can be characterized locally by their Gaussian
curvature. If $\kappa_{G}>0$ the surface locally bends in the same direction
at $P$, the geometry is that of a dome, Fig.3a. If $\kappa_{G}<0$, the surface
bends in opposite directions at the two principal directions. The surface
locally looks like a saddle, Fig.3b. If $\kappa_{G}=0$ (one of its principal
curvatures at least is 0) the surface is locally flat. Remarkably, this
characterization is purely intrinsic to the surface and does not depends on
its embedding in a higher dimensional space [9].

Figure 3: Principal curvatures on a surface When the vertical planes (in grey)
are rotated around the normal axis $\mathbf{n}$ (red arrow) with varying
angles $\theta$, the plane intersects the surface at curves $\gamma_{\theta}$.
The curvature vector $\mathbf{k}_{\theta}$ of these curves lies in the plane
made by the normal and the tangent vector, and is perpendicular to the tangent
vector $\mathbf{t}_{\theta}$. During rotation, the curvature intensity of
these intersection curves passes by a minimum and maximum value, which define
the principal curvatures. The corresponding directions are called the
principal curvature directions (indicated here by the black and dashed curves)
and are always perpendicular. (a) Surface with local positive Gaussian
curvature: the principal curvatures have the same sign. (b) Surface with local
negative Gaussian curvature: the principal curvatures have opposite signs.

###  3.2 Turtle state on a curved surface

On parametric surfaces, positions and directions can be thus specified
directly by defining positions and directions in the 2-D parameter space. For
instance, a position can be defined by providing a pair of coordinates
$P=(u^{1},u^{2})$ in the parametric space $\mathbb{U}$, while a direction at
this point may be defined by providing a vector $\mathbf{X}$ of coordinates
$X^{\alpha}$ in the local covariant basis, i.e.
$\mathbf{X}=X^{1}\mathbf{e}_{1}+X^{2}\mathbf{e}_{2}$. In $\mathbb{R}^{3}$, the
corresponding point coordinates are:

$(x^{1}(u^{1},u^{2}),x^{2}(u^{1},u^{2}),x^{3}(u^{1},u^{2})),$ (10)

where:

$x^{i}=J^{i}_{\alpha}u^{\alpha},$ (11)

where $J^{i}_{\alpha}$ are the components of the Jacobian $\mathbf{J}$:

$J^{i}_{\alpha}=\frac{\partial x^{i}}{\partial u^{\alpha}},$ (12)

that represent the best local linear approximation of $\phi$ in the
neighborhood of point $P$. This is called pushforward operator, and often
denoted $d\phi$. It is represented as a (3,2) matrix, made of the components
of vectors $\mathbf{e}_{1}$ and $\mathbf{e}_{2}$ arranged in two columns, and
maps vectors in the parameter space into corresponding vectors in the 3-D
space.

Because $\phi$ defines a bijective differentiable correspondence between the
points of $\mathbb{U}$ and those of the surface $\mathcal{S}$ that preserves
point neiborhoods, $x^{i}$ and $u^{\alpha}$ are often considered as the
coordinates of the surface point $P$ expressed in either $\mathbb{U}$ or
$\mathbb{R}^{3}$ respectively. Similarly vectors such as $\mathbf{X}$ and
$\mathbf{u}$ are mapped by the pushforward operator $d\phi$, and can be
interpreted as the ’same vector’ with coordinates expressed in two different
spaces, e.g. [7, p. 424].

Therefore the position and orientation of a turtle can be non-ambiguously
defined on the surface by specifying their coordinates in the parametric
space:

$\displaystyle P=$ $\displaystyle(u^{1},u^{2})$ $\displaystyle\mathbf{t}=$
$\displaystyle t^{1}\mathbf{e}_{1}+t^{2}\mathbf{e}_{2}.$

Together with the normal vector, the covariant basis form a local direct basis
of the 3D space at point $P$, $(\mathbf{e}_{1},\mathbf{e}_{2},\mathbf{n})$,
with basis vectors in the tangent plane or normal to it. As a consequence, one
can compute a unique $\\{\mathbf{H},\mathbf{L},\mathbf{U}\\}$ in
$\mathbb{R}^{3}$ on the curved surface. Assuming the turtle is at a position
$u^{\alpha}$ on the surface, and points in the direction of vector
$t^{\alpha}$ in the local tangent plane, the turtle’s head, $\mathbf{H}$, is
oriented in the direction of the tangent vector and is thus aligned along
$\mathbf{t}$. The up direction is chosen to be oriented along the surface’s
normal $\mathbf{n}$ and the turtle’s left arm points in the direction
perpendicular to both $\mathbf{U}$ and $\mathbf{H}$, Fig.4:

$\mathbf{H}=\displaystyle\frac{\mathbf{t}}{|\mathbf{t}|},\qquad\mathbf{U}=\mathbf{n},\qquad\mathbf{L}=\mathbf{U}\times\mathbf{H}$

Figure 4: Definition of the $\\{\mathbf{H},\mathbf{L},\mathbf{U}\\}$ frame on
a curved surface. (a) A direction $\mathbf{t}$ is defined at point $P$.
Together with the surface normal $\mathbf{n}$ at $P$, they define a local
reference frame that makes it possible to define (b) the turtle’s frame:
$\mathbf{H}$ is locally defined in the tangent plane aligned with the vector
$\mathbf{t}$ while the turtle’s upward direction $\mathbf{U}$ is imposed by
$\mathbf{n}$ and $\mathbf{L}$ is the direct vector product of $\mathbf{U}$ and
$\mathbf{H}$.

This make it possible to redefine the turtle state in a curved space. In
Euclidean space, this state was defined as:

$\displaystyle S=(x,y,z,H,L,U,\cdots,color,\cdots).$ (13)

In a curved surface, the turtle’s state now becomes:

$\displaystyle S=(\phi,u^{1},u^{2},t^{1},t^{2},\cdots,color,\cdots),$ (14)

where $\phi$ is the mapping
$\mathbb{U}\rightarrow\mathbb{R}^{3},x^{i}(u^{\alpha})$, defining the surface,
$(u^{1},u^{2})$ and $(t^{1},t^{2})$ are respectively the position and
orientation of the turtle in the parameter space.

#### SetSpace primitive

To construct a language based on the formalism of Riemannian L-systems we need
to associate language constructs with the main concepts introduced above. For
instance, to define the space within which the L-system will operate and the
turtle will move, the possibility to define parametric spaces must be
available in the language. The language primitive SetSpace will allow us to
select a specific parametric shape in a library of parametric surfaces
provided by the language and on which the movements of the turtle will take
place. This library contains standard geometric forms such as spheres, torus
or ellipsoids as well as more generic shapes such as surfaces of revolution,
sweeps, NURBS patches that make it possible to define more complex shapes.

⬇

1R = 2.

2Axiom: SetSpace(Sphere(R))

Once the parametric space is set, the turtle state takes the form indicated by
Equ.13. Then the turtle states are manipulated exactly as in the case of
classical Euclidean turtles, including the stacking and unstacking of turtle
states when reading L-strings, that push in and pop out states.

###  3.3 Moving straight in a curved space

Moving straight in a curved space means moving along geodesics. The definition
of geodesics relies on the notion of parallel transport, that specifies what
it means for vectors of a vector field to keep parallel as one moves along a
curve in the curved space.

#### Covariant derivative

To analyze the spatial variation of the vectors of a vector field as one moves
within a curved space, it is convenient to define a corresponding notion of
derivative. For this, the derivative of a vector at a point $P$ in a given
direction must not only integrate the change in coordinates of the vector, but
also the change of the local covariant basis. This leads to define a so-called
covariant derivative on vector fields.

To find an expression of covariant derivatives, let us consider the variation
of a vector in the tangent plane at a point $P$ along the coordinates lines in
the 3-D space:

$\begin{split}\frac{\partial\mathbf{X}}{\partial
u^{\alpha}}&=\frac{\partial(X^{\beta}\mathbf{e}_{\beta})}{\partial
u^{\alpha}}\\\ &=\frac{\partial X^{\beta}}{\partial
u^{\alpha}}e_{\beta}+X^{\beta}\frac{\partial e_{\beta}}{\partial
u^{\alpha}}.\end{split}$ (15)

This expression brings out the derivatives of the basis vectors, that can be
computed in the ambient space:

$\frac{\partial e_{\beta}}{\partial
u^{\alpha}}=\Gamma^{\gamma}_{\alpha\beta}e_{\gamma}+\Lambda_{\alpha\beta}\mathbf{n},$
(16)

where:

$\Gamma^{\beta}_{\alpha\gamma}=\frac{\partial\mathbf{e}_{\alpha}}{\partial
u^{\gamma}}.\mathbf{e}^{\beta},$ (17)

are the so-called Christoffel symbols (of second kind),
$\mathbf{e}^{\beta}=g^{\alpha\beta}\mathbf{e}_{\alpha}$ are the contravariant
basis vectors. The Christoffel symbols define the rate of change of the
covariant basis vectors in each direction of the local contravariant basis.
Similarly,

$\Lambda_{\alpha\beta}=\frac{\partial\mathbf{e}_{\alpha}}{\partial
u^{\beta}}.\mathbf{n},$ (18)

are the coefficients characterizing the variation of the covariant basis
vectors along the surface normal. Therefore, altogether we have:

$\frac{\partial\mathbf{X}}{\partial u^{\alpha}}=\left(\frac{\partial
X^{\gamma}}{\partial
u^{\alpha}}+X^{\beta}\Gamma^{\gamma}_{\alpha\beta}\right)e_{\gamma}+X^{\beta}\Lambda_{\alpha\beta}\mathbf{n}.$
(19)

Let us keep only the part of this expression lying in the tangent plane and
define:

$\begin{split}\nabla_{\alpha}\mathbf{X}&=\frac{\partial\mathbf{X}}{\partial
u^{\alpha}}-X^{\beta}\Lambda_{\alpha\beta}\mathbf{n}\\\ &=\left(\frac{\partial
X^{\gamma}}{\partial
u^{\alpha}}+X^{\beta}\Gamma^{\gamma}_{\alpha\beta}\right)e_{\gamma}.\end{split}$
(20)

This quantity is called the covariant derivative of $\mathbf{X}$ on the
surface. It corresponds to the orthogonal projection of the usual partial
derivative of the vector $\mathbf{X}$ in the Euclidean space onto the local
tangent plane, [66, p. 40]. It can be shown that $\nabla_{\alpha}\mathbf{X}$
has the properties of a derivative operator (linearity and product rule) and
that it defines an intrinsic differential operator on the surface, i.e. an
operator that depends only on quantities that can be measured within the
surface. Its components are denoted $\nabla_{\alpha}X^{\beta}$:

$\nabla_{\alpha}\mathbf{X}=(\nabla_{\alpha}X^{\beta})\mathbf{e}_{\beta}.$ (21)

More generally the covariant derivative of a vector $\mathbf{X}$ in the
direction of an arbitrary vector $\mathbf{Y}=Y^{\alpha}\mathbf{e}_{\alpha}$ is
defined as:

$\nabla_{\mathbf{Y}}\mathbf{X}=Y^{\alpha}\nabla_{\alpha}\mathbf{X}.$ (22)

This leads us to the definition of the notion of parallel transport. Let
$\gamma$ be a curve embedded within the surface defined from an interval
$\mathbb{I}\subset\mathbb{R}$ to the surface:

$\displaystyle\mathbb{I}$ $\displaystyle\rightarrow\mathcal{S}$ $\displaystyle
t$ $\displaystyle\rightarrow\gamma(t),$

and let $\mathbf{X}$ be a vector field on $\mathcal{S}$. $\mathbf{X}$ is
parallel transported along the curve $\gamma$ if at each point $\gamma(t)$ of
the curve:

$\nabla_{\mathbf{\dot{\gamma}}(t)}\mathbf{X}=0,$ (23)

where $\mathbf{\dot{\gamma}}(t)=\frac{d\gamma(t)}{dt}$ is the tangent to the
curve at point $\gamma(t)$. This means that the vector field keeps constant
seen from within the tangent planes when one moves in the direction of the
curve. More generally, two vectors fields $\mathbf{X}$ and $\mathbf{Y}$
parallel transported along a curve $\gamma(t)$ keep a constant scalar product
(i.e. a constant angle in the tangent planes along $\gamma$ as $t$ varies)
[7].

#### Geodesics

Consider a smooth curve $\gamma$ embedded in a curved space $\mathcal{S}$
(Fig.5a). We assume the curve $\gamma$ is parameterized by the arc-length
parameter $s$ using a smooth mapping $\gamma(s)=x^{i}(u^{\alpha}(s))$ from a
real interval $\mathbb{I}$ on the curved space $\mathcal{S}$, such that, as
$s$ varies, the point $\gamma(s)$ travels at a constant and unit velocity:

$\displaystyle\mathbf{t}=\frac{dP(s)}{ds}\quad\text{with}\quad|\mathbf{t}|=1.$
(24)

Figure 5: Curvature and geodesic. (a) Normal and geodesic curvature of a curve
lying on a surface. (b) Let $\mathcal{G}$ be the unique geodesic curve
starting at a point $P$ in the direction $\mathbf{X}$. $Q$ is the image on the
surface of $X$ by the exponential map at $P$: $Q=\exp_{P}(\mathbf{X})$.
Reciprocally, $\mathbf{X}=Log_{P}Q$.

The tangent vector being of constant norm, its rate of change along the curve,

$\mathbf{k}=\frac{d\mathbf{t}}{ds},$

reflects changes at every point $P$ in direction only and defines the
curvature vector $\mathbf{k}$ perpendicular to $\mathbf{t}$. The norm of this
vector,

$\kappa=|\mathbf{k}|,$

defines the curve’s curvature and indicates the intensity of the curve’s
bending in the 3D space at each point. Interestingly, on a curved surface
$\mathcal{S}$, the curvature vector $\mathbf{k}$ can be further decomposed
locally (Fig.5a). Let $\mathbf{n}$ and $T_{P}\mathcal{S}$ denote respectively
the normal and the tangent plane to $\mathcal{S}$ at a point $P$ (assuming a
local orientation of the surface). Then the curvature vector $\mathbf{k}$ can
be decomposed locally into a normal and a tangent components, $\mathbf{k_{n}}$
and $\mathbf{k_{g}}$,

$\mathbf{k}=\mathbf{k_{n}}+\mathbf{k_{g}}.$

$\mathbf{k_{n}}$ is the normal curvature (rate of change of the tangent vector
normal to the surface) and $\mathbf{k_{g}}$ is the geodesic curvature.
Geodesic curvature intensity indicates how the curve tends to bend locally in
plane, i.e. to deviate from a straight move in the tangent plane. Curves for
which the geodesic curvature is null at every point are called geodesics. As
for geodesics the change of tangent direction (acceleration) is purely normal
to the surface, a turtle living on the surface will not notice any lateral
movement when moving, and will have the impression to locally move straight.
Geodesics are the equivalent in curved spaces to straight lines in flat
euclidean spaces.

An important consequence of the definition is that, in the direction of the
tangent vector, the covariant derivative of the tangent vector of a geodesic
should have no components on the in-plane local basis vectors (otherwise there
would be a detectable geodesic curvature):

$\displaystyle\nabla_{\mathbf{t}(s)}\mathbf{t}(s)=\mathbf{k_{g}}=0$ (25)

This equation can be used to define geodesic curves on $\mathcal{S}$. Using
the definition of the tangent to the curve $\mathbf{t}$, Equ.24, expressed in
the covariant basis, i.e.
$\mathbf{t}=\frac{du^{\alpha}(s)}{ds}\mathbf{e}_{\alpha}$ and developing
Equ.25, one obtains a set of two second order, non-linear, coupled
differential equations, one for each value of $\alpha$ (e.g. [7, p. 106]):

$\displaystyle\frac{d^{2}u^{\alpha}}{ds^{2}}+\Gamma^{\alpha}_{\beta\gamma}\frac{du^{\beta}}{ds}\frac{du^{\gamma}}{ds}=0,$
(26)

These equations can be used to compute geodesic trajectories on the surface
from different perspectives, depending on the choice of boundary conditions
(see below).

##### Exponential maps.

Smooth surfaces have a remarkable property: At any point of the surface and in
a given direction, there exists a unique geodesic that passes through this
point and whose tangent at this point points in the given direction [73, p.
133]. This property makes is possible to define a map between vectors from the
tangent plane $T_{P}\mathcal{S}$ at $P$ and the points of the surface, called
the exponential map, e.g. [9, p. 287]. For this, let $P\in\mathcal{S}$ and
$\mathbf{X}\in T_{P}{\mathcal{S}}$ and denote $\mathcal{G}_{P}(t,\mathbf{X})$
the unique geodesic originating at $P$ with an initial velocity $\mathbf{X}$
(i.e. $\dot{\mathcal{G}}_{\mathbf{P}}(0,\mathbf{X})=\mathbf{X}$). Note that on
a geodesic parameterized by a parameter $t$, the norm of the velocity must
stay constant all along the curve, e.g. [66, p. 47]. The exponential map is
defined by [9, p. 288]:

$\displaystyle\exp_{P}{\mathbf{X}}=\mathcal{G}_{P}(1,\mathbf{X}),$

i.e. the exponential function returns the point reached after traveling on the
geodesic for a time unit, at constant speed $\|\mathbf{X}\|$. Therefore,

$\displaystyle Q=\exp_{P}{\mathbf{X}},$

is the unique point at a on the geodesic $\mathcal{G}_{P}(t,\mathbf{X})$ at a
unit time reach from $P$ when moving at constant velocity $\|\mathbf{X}\|$,
Fig.5b. At least in sufficiently small neighborhoods of $P$, this map is a
diffeomorphism [9, p. 288]. It is thus possible to define a reciprocal map,
called the logarithmic map at $P$, such that:

$\displaystyle\mathbf{X}=\log_{P}{Q},$

that returns, for any point $Q$ of the surface (in the region where $\exp_{P}$
is bijective), the direction $\mathbf{X}$ in the plane $T_{P}\mathcal{S}$
which initiates a geodesic from $P$ that passes through $Q$ and the geodesic
distance between $P$ and $Q$ has the norm of $\mathbf{X}$, e.g. [70, p.
17-30].

Exponential and logarithmic maps are essential tools in Riemannian L-systems
as they provide the natural concepts to simulate turtle movements locally,
i.e. in the local referential transported with the turtle during the movement.

#### Moving forward in a direction: an initial value problem

Consider a turtle positioned at a point $P$ on a curved surface, and heading
in direction $\mathbf{H}$. To advance the turtle by a distance $l$ on the
surface, one needs to compute the geodesic starting at $P=(u,v)$ in the
direction $\mathbf{H}=(p,q)$ over a length $l$ on the surface, Fig.6a. This
computation determines the new position of the turtle, corresponding to
$P^{\prime}=\exp_{P}(l\mathbf{H})$ and a new value of the head vector
$\mathbf{H}$ corresponding to the tangent of the geodesic at the destination
point $P^{\prime}$.

Figure 6: Different ways to specify straight displacement in a curved space.
(a) Forward algorithm: an initial value problem. (b) LineTo algorithm: a
boundary value problem

For this, we can integrate Equ.26 over a determined length $l$ to obtain the
points of the unique geodesic solving this initial value problem (IVP). A
classical strategy to solve an IVP system of second order differential
equations similar to Equ.26, consists of considering $\frac{du}{ds}=p$ and
$\frac{dv}{ds}=q$ as two new independent variables and rewrite Equ.26 as a
system of four coupled first order differential equations (see e.g. [22]):

$\displaystyle\begin{cases}\displaystyle\frac{du}{ds}=p\\\ \\\
\displaystyle\frac{dv}{ds}=q\\\ \\\
\displaystyle\frac{dp}{ds}+\Gamma^{1}_{11}p^{2}+2\Gamma^{1}_{12}pq+\Gamma^{1}_{22}q^{2}=0\\\
\\\
\displaystyle\frac{dq}{ds}+\Gamma^{2}_{11}p^{2}+2\Gamma^{2}_{12}pq+\Gamma^{2}_{22}q^{2}=0,\end{cases}$
(27)

with the initial conditions:

$P=\begin{bmatrix}u_{P}\\\
v_{P}\end{bmatrix},\mathbf{t}_{0}=\begin{bmatrix}p_{P}\\\ q_{P}\end{bmatrix}.$

The variables of Equ.26 have been explicitly rewritten with renaming
$u^{1}\rightarrow u$ and $u^{2}\rightarrow v$ to avoid over indexation in the
final Equ.27. This new system of first-order differential equations can be
integrated with classical ODE solvers. Note that for a given parametric
surface, the Christoffel symbols appearing in the equation are the only
quantities that need be computed. They can be evaluated by using Equ.17.

In a curved space, a forward step is specified by the primitive
$\texttt{F}(l)$ (move forward at the surface over a distance $l$) as in a
Euclidean space, except that the primitive will now make a call to a geodesic
solver that takes as an input the surface parametric equation, the current
position of the turtle $P=(u,v)$, the heading direction prescribed by the
current value of the turtle direction $\mathbf{H}$, both provided by their
coordinates in the parametric domain, as well as the number of elementary
steps $N$ requested to produce the line, Fig.6a. The solver returns a list of
$N$ points $P_{n}=(u_{n},v_{n})$, with their corresponding tangents
$\mathbf{t}_{n}=(p_{n},q_{n})$ on the geodesic that the turtle can use to move
forward, and draw a geodesic line on the surface as points of coordinates
$\\{P_{n}=x^{1}(u_{n},v_{n}),x^{2}(u_{n},v_{n}),,x^{3}(u_{n},v_{n})\\}_{n=0,\dots,N}$,
where $P_{N}=P^{\prime}=\exp_{P}(l\mathbf{H})$.

Altogether, in the language, all this integration process is hidden from the
user and results in very intuitive moves. The user can think of programming a
form as if moving in a Euclidean space. After the declaration of the
parametric surface to be used, moving forward by a distance $l$ length units
takes exactly the same form as moving forward with the turtle in a flat space:

⬇

1Axiom:

2 nproduce SetSpace(Sphere(1))

3 nproduce InitTurtle((0,0,1,0)) # Initializes turtle position at
u=0,v=0,p=1,q=0

4 nproduce F(0.5) # Draws a geodesic of length=0.5 from point (0,0) in
direction (1,0)

#### Moving forward to a given position on the curved space: a boundary value
problem

To prescribe turtle movements, an alternative method consists of specifying a
target point $P^{\prime}$ that the turtle should reach, from the current
position $P$, in the straightest possible way at its next step, Fig.6b. Here,
only a target point is prescribed instead of a direction and a length as in
the previous IVP case, and one seeks for the geodesics that corresponds to the
straightest path between the $P$ and $P^{\prime}$. This can be viewed as a
reciprocal problem, corresponding to evaluate the logarithmic map for a given
point $P^{\prime}$, i.e. find the direction $\mathbf{t}$ in $T_{P}\mathcal{S}$
and the corresponding geodesic that leads to $P^{\prime}$ from $P$, i.e. such
that $\mathbf{t}=\log_{P}P^{\prime}$.

Interestingly, the geodesic Equ.27, can also be used to solve this problem.
However, the problem is now constrained by the end-points, Fig.6b: given an
initial point $P$ and a target point $P^{\prime}$, find a geodesic that
connects these two points. For a smooth surface, there exists at least one
geodesic between two points. Finding such a geodesic belongs to the class of
boundary value problems (BVP). There are two main ways to solve these BVP
problems, using i) shooting methods or ii) improving progressively an initial
solution passing through the two endpoints, see e.g. [40].

##### Shooting method.

The first method can be formalized as an optimization problem by choosing a
shooting direction $\mathbf{t}$ as well as a length $l$, and integrating
Equ.27 over the length $l$ to find a geodesic and a length that would lead
exactly to the target point. Let $\gamma_{P,\mathbf{t}}$ be the unique
geodesic starting from $P$ in direction $\mathbf{t}$. For at least one
specific direction $\mathbf{t^{*}}$ ($\|{\mathbf{t^{*}}}\|=1$) this geodesic
passes through the target point $P^{\prime}$. Then for a certain optimal
shooting length $l^{*}$ we get:

$\exp_{P}(l^{*}\mathbf{t^{*}})=P^{\prime}.$

Then, by shooting in different directions and using different lengths, one can
evaluate the quality of each choice $(\mathbf{t},l)$ by computing a distance
$D$ between the point reached and the target $P$. The shooting problem can
thus be cast into an optimization problem:

$(\mathbf{t}^{*},l^{*})=\operatorname*{\arg\\!\min}_{\mathbf{t}\in
T_{P}(\mathcal{S}),l>0}D(\exp_{P}(l\mathbf{t}),P^{\prime}).$

The distance $D(.,.)$ should in principle be the distance on the surface. But
as computing this distance would already require the problem to be solved, in
practice, we use Euclidean distances in either the parameter space or the
embedding space that provide simple computations.

To solve this problem, classical numerical optimization techniques can be
used. Here we implemented this method using a nonlinear least-squares
algorithm with bounds on the variables [54]. The bounds on the variables make
it possible to fix the endpoints to the required values, $P$ and $P^{\prime}$
respectively.

The RiemannLineToShoot primitive can be called by passing the target point
$(u,v)$ coordinates as an argument:

⬇

1Axiom:

2 nproduce SetSpace(Sphere(1))

3 nproduce InitTurtle((0,0,0,1)) # Iniializes turtle position at
u=0,v=0,p=0,q=1

4 nproduce RiemannLineToShoot((1,0)) # Draw a geodesic between points (0,0)
and (1,0)

##### Geodesic residuals (GR) method.

In 1996, Maekawa [40] proposed an alternative optimization method to solve
this problem on parametric surfaces. Let $\Gamma(P,P^{\prime})$ be the set of
parametric curves with extremities fixed at $P$ and $P^{\prime}$ and $\gamma$
be a curve in $\Gamma(P,P^{\prime})$. The idea is to define a quantitative
criterion $C(\gamma)$ to assess how much $\gamma$ departs from a geodesic,
($C(\gamma)=0$ means the curve is a geodesic) and then to minimize this
criterion:

$\gamma^{*}=\operatorname*{\arg\\!\min}_{\gamma\in\Gamma(P,P^{\prime})}C(\gamma),$

where $\gamma^{*}$ is the sought geodesic between end-point $P$ and
$P^{\prime}$. To define the criterion, let us remark that Equ.27 can be
written in the form:

$\frac{d\mathbf{Q}}{ds}=G(\mathbf{Q},s),$ (28)

where $\mathbf{Q}$ is the vector $[u,v,p,q]^{T}$ corresponding to the
concatenation of the position and velocity of a point on the geodesic. In
addition, the geodesic curve must respect the boundary conditions:

$\mathbf{Q}_{initial}=[u=u_{P},v=v_{P},-,-]^{T}\quad\textrm{and}\quad\mathbf{Q}_{final}=[u=u_{P^{\prime}},v=v_{P^{\prime}},-,-]^{T},$
(29)

where $(u_{P},v_{P})$ and $(u_{P^{\prime}},v_{P^{\prime}})$ are the
coordinates of the two end-points and a dash $-$ means that these values are
unconstrained at the boundary. This equation is a first order differential
equation that suggests that recurrence relations exist binding the variables
along the curve as $s$ varies. This can be made explicit by discretizing the
curve into $m-1$ segments bounded by points of curvilinear abscissa
$[s_{k},s_{k+1}],k=0,\dots,m-1$. Denoting $\mathbf{Q}_{k}=\mathbf{Q}(s_{k})$
and $\mathbf{G}_{k}=\mathbf{G}(\mathbf{Q}_{k},s_{k})$, and using the
trapezoidal rule to approximate the derivative, one obtains the following
recurrence relations:

$\frac{\mathbf{Q}_{k}-\mathbf{Q}_{k-1}}{s_{k}-s_{k-1}}=\frac{1}{2}(G_{k}+G_{k-1}),\qquad\forall
k=1,\dots,m-1,$

with the boundary conditions defined in Equ.29 becoming
$\mathbf{Q}_{0}=\mathbf{Q}_{initial}$ and $\mathbf{Q}_{m}=\mathbf{Q}_{final}$.
Note that the distance $s_{k}-s_{k-1}$ between curvilinear abscissa needs to
be approximated by a Euclidean distance in $\mathbb{R}^{3}$. However, provided
that the segments are sufficiently small compared to the local curvature, this
approximation is in general accurate. From this equations, let us define the
residuals:

$R_{k}=\frac{\mathbf{Q}_{k}-\mathbf{Q}_{k-1}}{s_{k}-s_{k-1}}-\frac{1}{2}(G_{k}+G_{k-1}),\qquad\forall
k=1,\dots,m-1,$ (30)

and

$R_{0}=[u_{0}-u_{P},v_{0}-v_{P}]^{T}\qquad\textrm{and}\qquad
R_{m}=[u_{m}-u_{P^{\prime}},v_{m}-v_{P^{\prime}}]^{T}.$ (31)

For points falling exactly on a geodesic, these residuals should be 0. The
problem can thus be turned into an optimization problem where we seek for the
$\mathbf{Q}_{k}$ and $s_{k}$ that cancel the residuals:

$\displaystyle\begin{split}R_{k}=&[0,0,0,0]^{T}\qquad k=1,\dots,m-1,\\\
R_{0}=&R_{m}=[0,0]^{T}.\end{split}$ (32)

This defines a non-linear system of $4m=4(m-1)+2+2$ equations for
$4m=4(m-1)+2+2$ variables. Here again, such a system can be solved using
classical methods to find the roots of a function. Following Maekawa [40] who
gives an explicit expression for the Jacobian of the residuals, we used a
Newton method that remains very efficient if a good initial guess of the
solution can be found. In our case, a natural initial guess corresponds to the
linear interpolation in the parameter space of the endpoint variables
$(u_{P},v_{P})$ and $(u_{P^{\prime}},v_{P^{\prime}})$.

In Riemannian L-systems, the GR method can be called using a RiemannLineTo
primitive and passing the target point $(u,v)$ coordinates as an argument.
This computes the geodesic between the current turtle position and the target
point argument, and moves the turtle’s position at the target point, with a
heading direction $\mathbf{H}$ tangent to the geodesic extremity.

⬇

1Axiom:

2 nproduce SetSpace(Sphere(1))

3 nproduce InitTurtle((0,0,0,1)) # Iniializes turtle position at
u=0,v=0,p=0,q=1

4 nproduce RiemannLineTo((1,0)) # Draw a geodesic between points (0,0) and
(1,0)

#### Geodesic distances on the surface

Distances can be defined between points of a curved space by using the length
of the shortest path between two points. In addition to being the straightest
lines between points, geodesic also have the property of being local minima of
the length of trajectories between any two points. We can thus define a notion
of distance between points of the surface by integrating the length of the
different segments composing a geodesic. Let us call $D_{\mathcal{S}}$ this
geodesic-based distance on the surface $\mathcal{S}$ and call
$\\{(u_{k},v_{k})\\}_{k=0,m}$ the points in the parameter space defining a
geodesic between $A$ and $B$, then:

$D_{\mathcal{S}}(A,B)=\sum_{k=1}^{m}\|\mathbf{x}(u_{k},v_{k})-\mathbf{x}(u_{k-1},v_{k-1})\|,$
(33)

where $\|.\|$ is the usual $L_{2}$ norm in $\mathbb{R}^{3}$. Here, we assume
that the distance in $\mathbb{R}^{3}$ approximates sufficiently well the
geodesic length of short segments on the surface. Note that in general the
geodesic between 2 points of the surface is not necessarily unique. In
principle, if several geodesic exist one should base the above distance on the
one that minimizes their lengths.

The distance between two points can be computed by using the function
geodesic_distance_to_point, as soon as the turtle state is available (mainly
in interpretation rules).

⬇

1 dist,_,_ = geodesic_distance_to_point(turtle.space, (u,v),(ut,vt))

###  3.4 Turning on curved surfaces and holonomy

Riemannian spaces and in particular parametric surfaces, are by definition
smooth spaces that can be locally approximated by Euclidean spaces. As turning
(i.e. rotating) in a Riemannian space is a purely local operation, i.e. an
operation that takes place in the tangent plane associated with the current
position, it is thus no wonder that turning in such curved spaces is
essentially equivalent to turning in Euclidean spaces. However, one must be
cautious as rotation transformations must usually be operated in orthonormal
bases. This is not the case in our Riemannian spaces as the local basis, e.g.
the covariant basis, is not in general orthonormal. Before applying a
rotation, one must thus move temporarily to some local orthonormal basis and
move back to the original basis afterwards.

Figure 7: Turning on a surface. (a) On a geodesic trajectory (black curve),
the turtle is instructed to turn by an angle $\theta$ at a point $P$. As the
tangent $\mathbf{t}$ (black) before the turn is expressed in the local
covariant basis (red), a orthonormal basis (not shown) must be computed to
perform the rotation, leading to a new tangent vector $\mathbf{t^{\prime}}$
expressed in the orthonormal basis. Then, the new tangent vector is expressed
in the covariant basis and moves along geodesics can continue. (b) in this
way, the turtle can draw curved polylines on the surface, by alternating
geodesic segments (in black) and rotations (angles $\theta$ in blue).

The basis vectors $\mathbf{e^{\prime}}_{\alpha}$ of the new coordinate system
can be related to the old basis vectors $\mathbf{e}_{\beta}$ by a matrix
$\mathbf{G}=\\{G_{\alpha}^{\beta}\\}$ such that,
$\mathbf{e^{\prime}}_{\alpha}=G_{\alpha}^{\beta}\mathbf{e}_{\beta}$ and
$\mathbf{e^{\prime}}_{\alpha}.\mathbf{e^{\prime}}_{\beta}=\delta_{\alpha\beta}$.
Then to rotate in the local tangent plane a vector $\mathbf{t}$ with
components $[t^{1},t^{2}]^{T}$ in the covariant basis by an angle $\theta$,
one must first transport the components of $\mathbf{t}$ in a local orthonormal
basis, then do the rotation yielding the vector $\mathbf{t^{\prime}}$ with
components expressed in the orthonormal basis, and finally get back to the
original, non-orthonormal, basis where $\mathbf{t^{\prime}}$ has components
$[t^{\prime 1},t^{\prime 2}]^{T}$, Fig.7. If $\mathbf{R}_{\theta}$ denotes the
rotation matrix by an angle $\theta$ in the orthonormal basis
$\mathbf{e^{\prime}}_{\alpha}$, then:

$\begin{bmatrix}t^{\prime 1}\\\ t^{\prime
2}\end{bmatrix}=\mathbf{GR_{\theta}G^{-1}}\begin{bmatrix}t^{1}\\\
t^{2}\end{bmatrix}.$ (34)

Hence for turning by an angle $\theta$, one must first compute locally a
matrix $\mathbf{G}$ (here we used a classical Gram-Schmidt orthogonalization
process [16].) and then rotate turtle’s head $\mathbf{H}$ and left arm
$\mathbf{L}$ according to Equ.34, while keeping the $\mathbf{U}$ vector
unchanged in the direction of the surface normal. All these operations are
hidden in the language and the user only specifies rotations by providing
angles (expressed in degrees):

⬇

1Axiom:

2 nproduce SetSpace(Sphere(1))

3 nproduce InitTurtle((0,0,0,1))

4 nproduce F(1) +(30) F(1) +(45) F(2)

##### Parallel transport

The possibility to turn while moving in a curved space is essential to reveal
the curved nature of the space in which the movement takes place, e.g. [3, p.
301] and [1]. To see this, let us use the notion of parallel transport
introduced above and see how it can be illustrated using the turtle movements.
A vector is parallel transported along a geodesic, its orientation with
respect to the geodesic tangent and its norm does not change as one moves
along the geodesic.

Figure 8: Holonomy and parallel transport using Riemannian Lsystems. (a-b)
parallel transport of a vector along a polygon made of geodesics. The vector
is initially tangent to the first geodesic, then is perpendicular to the
tangent on the second geodesic, then points backward on the third geodesic
segment. (d-c) parallel transport of a vector not tangent to the first
geodesic. The vector keeps a constant angle with the tangent vector, but this
angle changes each time the turtle turns. (e-f) parallel transport along a
curve that is not a geodesic: the angle between the transporting curve and the
transported vector varies continuously.

Consider the classical example of a walker moving on a sphere and transporting
a vector. The walker starts at the equator in the direction of the North pole
and holds her vector pointing ahead (e.g. [1], Fig.8a-b red arrow pointing
upward at green point). While walking straight, the walker keeps the vector
always positioned identically with respect to her. She is parallel-
transporting the vector with her. The vector is thus aligned with the tangent
of the curve and keeps aligned all the way through. At the pole, the vector is
thus horizontal and points to the back. The walker then turns to the left
without turning the vector. She then continues her path while still parallel
transporting the vector, Fig.8a. As the orientation of the vector with respect
to the geodesic tangent does not change along the way, the vector arrives at
the equator pointing to the back. The walker then turns again by 90 degrees
counter clockwise, still without moving the transported vector orientation,
and continues her trip along the equator back to its starting point. When
arriving, the transported vector now points west in the initial tangent plane.
Still without changing the orientation of her transported vector, the walker
makes a final turn by 90 degrees to get back to its exact starting orientation
and completely close the loop. One can observe that, although during the trip,
orientation of the transported vector never got modified, its final
orientation (to the west) does not match its original orientation (to the
north).

This property of implicitly rotating parallel transported vectors along closed
trajectories is characteristic of curved spaces. It is called holonomy. The
difference of angle between the vectors before and after the parallel
transport along a loop is called the angle defect. In flat spaces, the angle
defect is $0$ for any loop and any parallel-transported vector. In curved
spaces, the angle defect is tightly linked with the curvature of the space,
e.g. [1].

Holonomy and angle defect can be vizualized with Riemannian L-system. For this
the turtle keeps track of the accumulated rotation angle, $\alpha$, since an
origin point on a turtle path. The origin point is defined using the module
ParallelTransportReset that reiniliatizes the cumulated turtle rotation at any
moment ($\alpha=0$), and thus defines a new origin for parallel transport at
the current turtle’s position, say $(u_{0}=0,v_{0}=0)$. For a unit vector
$\mathbf{v}$ positioned in the tangent plane at the origin point and making an
angle $\beta$ with the turtle’s head at the origin, we can thus compute and
draw its corresponding transported vector at any subsequent position of the
turtle. This is done by using the module
ParallelTransportedArrow(vect_angle,vect_size), with $vect\\_angle=\beta$, and
$vect\\_size$ being a scaling factor for vizualization. At the new position,
the parallel transported vector $\mathbf{v}$ makes an angle $\beta^{\prime}$
with the current turtle’s head in the current tangent plane corresponding to
the original angle corrected by the cumulated angle (i.e. $\alpha$) by which
the turtle turned on the path to the current point:

$\beta^{\prime}=\beta-\alpha$

The module ParallelTransportedArrow(vect_angle,vect_size) thus draws the
transported vector $\mathbf{v}$ as a vector making an angle $\beta^{\prime}$
with the turtle’s head at the current turtle’s position.

To illustrate how to program this with Riemannian L-systems, let us start with
a program tracing a triangle at the surface of a sphere as in Fig.8c:

⬇

1Axiom:

2 nproduce SetSpace(Sphere(1))

3 nproduce InitTurtle((0,0,0,1))

4 nproduce F(R*pi/2) +(90) F(R*pi/2) +(90) F(R*pi/2) +(90)

This program can be modified to show the different stages of a transported
vector along the triangular path. For this, after setting the origin point of
parallel transport at 0 using ParallelTransportReset, we insert the module
ParallelTransportedArrow() at each vertex positions (Fig.8c-d):

⬇

1alpha0 = -30 # in degrees

2R = 1

3Axiom:

4 nproduce SetSpace(Sphere(R))

5 nproduce InitTurtle((0,0,0,1))

6 nproduce ParallelTransportReset

7 nproduce ParallelTransportedArrow(alpha0,0.5)

8 nproduce F(R*pi/2)

9 nproduce ParallelTransportedArrow(alpha0,0.5)

10 nproduce +(90) F(R*pi/2)

11 nproduce ParallelTransportedArrow(alpha0,0.5)

12 nproduce +(90) F(R*pi/2)

13 nproduce ParallelTransportedArrow(alpha0,0.5)

Listing 5: Parrallel transport (see Fig.8.c-d)

Fig.8e-f, illustrating the transport of a vector along a non-geodesic curve
can be obtained in a similar manner by making short segments of size $dl\ll
R$, each followed by small rotations $d\alpha\ll\pi$ and drawing the
transported vector at each incremental step.

##### Interpretation of turtle movements in terms of differential operators.

According to the definition of parallel transport, we see that the head vector
of the turtle $\mathbf{H}$ is always parallel transported by a F statement as
it remains parallel to the tangent vector of the geodesic trajectory produced
by F. Hence, parallel transport of a vector $\mathbf{X}$ in the current
tangent plane $T_{P}\mathcal{S}$ of the turtle, through a F statement, is the
vector $\mathbf{X^{\prime}}$ such that, if $\alpha$ denotes the angle between
$\mathbf{X}$ and $\mathbf{H}$ in $T_{P}\mathcal{S}$, $P^{\prime}$ and
$\mathbf{H^{\prime}}$ are respectively the new position and head of the turtle
after the execution of the F statement, then $\mathbf{X^{\prime}}$ is the
vector in $T_{P^{\prime}}\mathcal{S}$ that makes an angle $\alpha$ with
$\mathbf{H^{\prime}}$ and such that $\|\mathbf{X^{\prime}}\|=\|\mathbf{X}\|$.

Combined together, the Riemanian definition of a forward (F) and rotate (+)
statement, implement the exponential map at the current turtle’s position $P$.
Let us take the turtle’s orientation $\mathbf{H}$ as a reference orientation
in current turtle’s tangent plane. Then any vector $\mathbf{X}$ in this plane
can be defined by a rotation $\alpha$ with respect to $\mathbf{H}$, and a
scaling factor $\|\mathbf{X}\|=l$. Then the instruction,

⬇

1 nproduce +(alpha) F(l)

moves the turtle at the point $P^{\prime}=\exp_{P}(\mathbf{X})$ with a new
direction $\mathbf{H}^{\prime}$ at $P^{\prime}$ that is parallel transported
from the old direction $\mathbf{H}$ at $P$ along the geodesic from $P$ to
$P^{\prime}$.

Reciprocally, given the current turtle’s position $P$ and a different point
$P^{\prime}$ on the surface, the definition of the statement RiemannLineTo
implements the logarithmic map at the turtle’s position. The instruction:

⬇

1 nproduce RiemannLineTo(P_prime)

computes a geodesic as a ordered list of points on the geodesics
$\\{(u_{n},v_{n},p_{n},q_{n})\\}_{n=0,\dots,N}$, where
$P_{N}=P^{\prime}=\exp_{P}(l\mathbf{H})$. The vector
$\mathbf{X_{0}}=(p_{0},q_{0})$ belongs to the tangent plane at $P$ and
represents the initial direction of the velocity in the parameter space, and
thus $\mathbf{X_{0}}=\log_{P}(P^{\prime})$.

##### Closed polygons on a surface

As illustrated by holonomy, various usual geometric properties valid in
Euclidean spaces are no longer valid in curved spaces. In particular, the sum
of the inner angles of a triangle is not 180 degrees nor even a constant
number in general. This makes it difficult to draw closed polygonal curves
just based of local operations such as moving forward by a length $l$ or
turning by a certain angle $\alpha$, Fig.9a. In general, is not possible to
map flat space onto a curved space map while preserving both the angles and
lengths. This means that we have to accept to drop some properties of polygons
when mapping them in a curved space. However, it can be practically important
in some cases to draw close figures on the surface and even to fill them with
a particular color or texture.

This issue related to holonomy can be circumvented in different ways in
Riemaniann L-systems. A first option is to draw a polygonal line using the F
module (thus making geodesic segments) Fig.9a, and to close the polygon by
creating a final geodesic segment using a RiemannLineTo module from the
current position of the turtle to its initial position. However, this solution
may induce undesired biases on the length and orientation of the final segment
if the total curvature enclosed by the polygonal line is significant.

A second option consists of tracing the polygonal figure in the parametric
space rather than on the surface directly, and then project all the points of
this polygonal line on the surface Fig.9b. This can be considered as an
indirect interpretation of the turtle instructions, which are momentarily
executed in the parameter space rather than on the surface. This can be
carried out by using the modules StartIndirectInterpretation and
StopIndirectInterpretation.

The backside of this indirect interpretation is that the segments are in
general not geodesics. It is possible to define a spherical square as a
polygon made of geodesic edges with equal lengths and equal inner angles
between consecutive edges (although different from 90 degrees). In Riemannian
L-systems one would position the turtle at the center of the square, then send
geodesic rays of equal length $d$ from this point with an angle of 90 degrees
between each other (Fig.9c thin blue lines). This defines four points at the
same distance of the square center, which can then be joined by a geodesic
using a RiemannianLineto primitive. This procedure can be used on more chaotic
surfaces (Fig.9d), but with an additional loss of geometric symmetries of the
square.

Figure 9: Drawing closed polygons on curved surfaces. (a) failure to close a
square by forcing consecutive sides to be at 90 degrees from one another on a
curved surface. (b) Square drawn in the parameter space, and then pushed on
the surface (c) Alternative geometric construction of the square preserving
the right angle at the intersection of the diagonals and their length. (d) the
construction in c can be used on more chaotic surfaces.

##### Drawing smooth curves on a curved surface

In the Euclidean space, important families of parametric curves can be
constructed by the use of straight lines. This is the case for instance of
basis spline curves, or B-Splines, whose points are a weighted linear
combinations of so called control points, [53]. This linear combination is
based on piece-wise polynomial basis functions whose degree controls the
smoothness of the curve. An approximation of such curve can be achieved
through algorithms that recursively subdivide the control polygon (polygon
formed by the control points), such as the de Casteljau [53] or the Lane-
Riesenfeld [35] algorithms. Here, we show how such algorithms can easily be
extended with Riemannian L-systems to control smooth curves in curved spaces
(Fig. 10.a-h).

Starting from the initial control polygon made up of geodesics segments, the
number of controls points is first doubled by inserting new control points in
the middle of each segment (duplication operation) (Fig. 10.b). The centroid
of the geodesic is estimated as the point of the geodesic at equal geodesic
distance between the two control points. Then, each initial control point is
moved toward the midpoint of the adjacent segment (move operation)(Fig.
10.c-d). Each subdivision step is thus composed of duplication and move
operations. Subdivision steps are reiterated until the desired level of
approximation of the B-Spline curve is reached (Fig. 10.f-h).

Let us illustrate how this principle can be used to generate an oak leaf-like
shape using Riemannian L-systems. First, a set of control points is specified
as the terminal nodes of a simple tree structure using the BSplinePoint
module, encapsulated within StartBSpline and EndBSpline modules 6. These
control points must be provided in an order that respects the curve
parameterization. For this, the left branches are first generated (line 17),
followed by the primary branch (line 18) and subsequently the right ones (line
19). The subdivision algorithm encoded in the BSpline primitive of the system
(line 6,8) then generates the outline of an oak-like leaf. Depending on the
embedding space, the shape is drawn on a flat surface 10i, an ellipsoid
Fig.10j or a bumpy ellipsoid respectively in Fig.10k. The representation of
the leaf contour dynamically adjusts to the local irregularity of the
underlying space. In particular, on the bumpy ellipsoid, the main axis of the
skeleton is deflected toward the right, due to the specific curvature of the
space.

Various types of shapes can be drawn on curved spaces using the same principle
as illustrated by Fig.10l-n. Here, the control polygon represents a salamander
sketched in Escher style.The original control polygon (Fig. 10l) is
recursively subdivided to achieve approximation of B-Spline curves of degree 2
and 8 respectively (Fig. 10m-n). The higher the degree, the smoother the curve
and more rough the details.

Figure 10: B-splines on a curved surface. (a) Initial quadrilateral control
polygon on a sphere. (b) Duplication operation with new control points (in
purple) inserted in the midpoint of each segment. (c-d) All original control
points (in red) are moved toward the midpoint of their adjacent segments (in
purple). For a B-Spline of degree 3, 2 move operations are applied. (e)
Resulting control polygon after a complete subdivision step. (f-g-h)
Successive control polygons (in green) after 1, 2 and 3 subdivision steps
respectively. (i-k) The B-Spline curve is defined by control points positioned
using the Riemannian L-system 6 that generates a simple branching structure.
Resulting curve on (i) a flat surface, (j) an ellipsoid and (k) a bumpy
ellipsoid. (l) The control polygon of a salamander shape. (m-n) The
interpretation of the control polygon as B-Splines of degree 2 and 8
respectively. The degree of the curve controls the number of control points
that influence each point of the curve.

⬇

1maxlength = 12

2dl = 1.5

3

4Axiom:

5 nproduce SetSpace(EllipsoidOfRevolution(5,8)) InitTurtle([0,-1,0,1])

6 nproduce StartBSpline(2) # The degree of the B-Spline is given as parameter

7 nproduce BSplinePoint() [ A(0) ] BSplinePoint()

8 nproduce EndBSpline()

9

10derivation length: 8

11production:

12

13A(length) :

14 if length < maxlength:

15 length += dl

16 lateral_length = maxlength*lateralratios(length/maxlength)

17 nproduce F(dl)

18 nproduce [+(60) nproduce F(0.1+lateral_length) BSplinePoint()]

19 nproduce [ F(0.1) A(length)]

20 nproduce [-(60) F(0.1+lateral_length) BSplinePoint()]

21

22interpretation:

23A(l) \--> BSplinePoint()

Listing 6: B-Spline curve built by positionning B-Spline control points at the
end of a simple tree structure (see Fig.10j). The lengths of the lateral
branches depend on a graphically defined function called lateralratios.

##  4 Freely growing forms on curved surfaces

Using the turtle primitives introduced above to move on curved surfaces, we
now explore how to program the development of filamentous or branching forms
on curved surfaces using Riemannian L-systems.

###  4.1 Geodesic trajectories

In the absence of any reason for deviating from a straight movement, a mobile
moving in a curved space would naturally follow geodesics. On the sphere for
example, starting at a point $P$ and heading in a direction $\mathbf{v}$ will
produce a geodesic that corresponds to the great circle passing through $P$ in
the direction $\mathbf{v}$. In more complex spaces, geodesics may show
remarkable properties that can lead to non-intuitive patterns.

Let us consider for example geodesic trajectories on an ellipsoid of
revolution, starting at equatorial points in directions given by different
initial angles above the equator $\alpha$, Fig.11a. In Riemannian L-systems
these geodesics can be simulated with the simple program below:

⬇

1alpha = -45 # in degrees

2ra = 1. # dimensions of the ellipsoid

3rb = 0.5

4lg = 10 * (2*pi*ra) # length of the geodesic

5Axiom:

6 nproduce SetSpace(EllipsoidOfRevolution(ra,rb))

7 nproduce InitTurtle((0,0,1,0)) # head points in the horizontal direction
[1,0]

8 nproduce +(alpha) # initial inclination with respect to horizontal

9 nproduce F(lg) # draws a geodesic

Listing 7: Geodesic on an ellipsoid of revolution

Here, only an axiom is defined, with no production rules, and directly
produces modules that can be interpreted by a Riemannian turtle (no need of
interpretation rules either). The axiom first triggers the use of Riemannian
L-systems by defining a curved space in which the turtle should operate (line
6), then its initial position and reference orientation (pointing right in the
horizontal direction) in this space (line 7), and finally after a local
rotation with respect to the reference orientation (line 8) a geodesic line of
length lg is drawn from the initial position heading at angle alpha below the
horizontal direction (line 9). By varying alpha, we can observe the behavior
of geodesics on an ellipsoid, Fig.11a. They tend to stay in an equatorial band
whose width depends on the inclination of the initial direction with respect
to the equator: the smaller the initial inclination $\alpha$, the narrower the
band. Interestingly, depending on the initial inclination angle, the geodesic
will more of less densely cover this equatorial band.

Figure 11: Geodesics on surfaces of revolution. (a) Geodesic on an ellipsoid
of revolution, with a length 10 x equator circumference (equator indicated in
blue), and with different initial orientation (green arrow): from top-left to
bottom-right: $60,45,44,43.3,43.2,43,30,10$ degrees inclination with respect
to equator. (b) Comparison of the behavior of close geodesic trajectories in
spaces with positive (left: sphere) and negative (right: pseudo-sphere)
curvatures. In both examples, geodesics start with parallel orientation (red
arrow). (c) Geodesics in a space with negative Gaussian curvature starting
with varied initial orientations. Geodesics are all the more deflected by the
space curvature that they pass closer to the center of the shape.

The curvature of the space can in general be characterized by looking at how
neighboring geodesics behave, e.g. [7, p.144]. Geodesics initially parallel
will tend to converge in spaces with positive Gaussian curvature and diverge
in space with negative Gaussian curvature, Fig.11b. It can be shown that this
property, called geodesic deviation, characterizes locally space curvature:
the rate at which parallel neighboring geodesics get closer or farther away
from each other is proportional to the local curvature [7].

Fig.11c shows an example of geodesics on a surface with negative curvature.
Several geodesics displayed with different colors are emitted in slightly
different directions toward the center of the surface (central hole). The
closer they pass to the axis of the central hole, the more the geodesics are
deflected. Geodesics passing too close to the hole are even reflected back
(orange curve).

These examples are actually the consequences of a these behaviors are the
consequences of a mathematical theorem about geodesics on surfaces of
revolution, known as the Clairaut’s theorem [9, p. 259] for Clairaut
parameterizations of general surfaces [50, p. 353]):

######  Theorem (Clairaut).

Consider a geodesic on a surface of revolution and denote $\alpha$ the
inclination angle of the a geodesic at a point $P$ with respect to the
latitude passing by $P$ and $r$ the radius of revolution at $P$. Then, along
the geodesic we have:

$r\cos\alpha=c_{0},$ (35)

where $c_{0}$ is a real constant.

Therefore, choosing an initial point and direction for a geodesic determines a
value $c_{0}$ that will remain subsequently constant along the geodesic
trajectory. As one moves along the geodesic, $r$ varies, and $\alpha$ must
vary accordingly to respect Equ.35. If $r$ decreases, $\cos\alpha$ must
increase. At some point, if $r$ continues to decrease, $\cos\alpha$ reaches
the value 1 (i.e. the curve is tangent to the local latitude at $P$), and the
trajectory will be reflected back to increasing values of $r$ to avoid $r$
decreases more and to keep up with Equ.35. This means that the geodesics with
this $c_{0}$ will be trapped in a part of the surface of revolution. For the
ellipsoid (Fig.11a) for example, the geodesics keep in an equatorial band
whose width is determined by the value of $c_{0}$ at the initial point and for
the chosen initial direction. For the surface with negative curvature
(Fig.11c), trajectories passing close to the central hole have a radius which
markedly decreases and thus induce a reversal of the angle $\alpha$ variation
at some critical radius, thus bouncing back the geodesic curve.

These examples illustrate how Riemannian L-systems can help explore
mathematical properties in differential geometry with the efficiency and
simplicity provided by a high-level programming language (see List.7).

###  4.2 Turning and branching

##### Fractals.

Being particular cases of manifolds, smooth parametric surfaces are locally
flat and look like a Euclidean plane. However, at larger scales, the surface
curvature is not negligible and can be revealed by trajectory holonomy.
Fractals living on curved surfaces exhibit various levels of detail all along
their entire structure. We may expect that the fine details of the structure,
typically much smaller than the local surface principal radii of curvature,
are not affected by the surface curvature. However, at coarser scales (of the
order of magnitude of the radii of curvature), the fractal form must be
affected in proportion of the local curvature. This is illustrated by the
series of prefractal forms in Fig.12. The reference prefractal forms converge
toward the well known von Koch flake in a flat space Fig.12a. When traced on a
sphere, we see that the curve does not close any more due to holonomy at large
scale that tend to fold the curve faster than in the flat reference case,
12.b. However, one can notice that the curve pattern is almost not affected at
fine scales where one easily locally recognizes the pattern of the von Koch
curve. Interestingly, if the radius of the sphere is reduced (the curvature is
increased), the overall topology of the flake is dramatically affected, and
contains only four main arms Fig.12c, instead of six in the reference case.
Similarly, Fig.12d, shows a von Koch curve growing on a torus. when the curve
size is small compared to the two main radii of curvature of the torus (left),
the curve form is not affected significantly. However, as the curve grows, the
overall curve topology is markedly deformed (right).

These figures have been obtained with the program listed below (List.8). One
can observe that the modification of the code with respect to the reference
flat case (List.3) is minimal: only the specification of the turtle space has
been added in the axiom. The rest of the code is unchanged. For the torus, one
just need to change line 3 by nproduce SetSpace(Torus(1,0.3)).

Figure 12: Using the turtle to draw fractals on curved spaces. (a) prefractal
sequence of the von Koch Curve in a flat space (b) Prefractal sequence
obtained by the same procedure as in a on a sphere of radius 1, and (c) on
sphere of radius 1/2. (d) von Koch curves with increasing step length on a
torus.

⬇

1Axiom:

2 nproduce SetSpace(Sphere(radius=1))

3 nproduce F(1)-(120)F(1)-(120)F(1)

4derivation length: 5

5production:

6F(x) : nproduce F(x/3.0)+F(x/3.0)--F(x/3.0)+F(x/3.0)

Listing 8: Fractal curve (von Koch flake) on curved surfaces (see Fig.12b-d)

##### Branching patterns.

Like in classical L-systems, branching patterns can be created easily with
Riemannian L-systems using well-formed strings of modules. The turtle’s
interpretation mechanism is entirely conserved but runs on Riemannian states
(Equ.14) instead of the classical turtle’s states used in Euclidean geometry
(Equ.13). The turtle thus draws branching patterns whose segments are
geodesics, and branching angles are defined either explicitly by the
programmer in the local tangent plane at the bifurcation point using a +
statement before an F, Fig.13a-c or implicitly inferred through a
RiemannianLineTo primitive. The generic program to produce these trees is as
follows (List.9):

⬇

1N = 7 # Depth of the tree recursion

2iangle = 45 # Insertion angle

3ilen = 0.2 # Length of a segment between two branches

4

5Axiom:

6 nproduce SetSpace(Sphere(radius=1))

7 nproduce InitTurtle([0.,0.,0.,1.]) # turtle’s head pointing upwards

8 nproduce A(0)

9

10derivation length: N

11production:

12A(n) :

13 a = iangle if n % 2 else -iangle

14 if n<N:

15 nproduce F(ilen) [+(a)A(n+1)] A(n+1)

Listing 9: Tree patterns on surfaces (see Fig.13a-c)

As before, to change the space (Sphere, Pseudo-sphere, …) and the initial
position of the turtle in the space, only lines 6 and 8 have to be updated.
The rest of the code remains unchanged.

Branching patterns, that are made of geodesic segments, are affected
differently by the surface depending on both its extrinsic and intrinsic
geometric properties, Fig.13a-c. Branches tend to get more dense on the sphere
(constant positive Gaussian curvature) Fig.13a, less dense on the pseudo-
sphere (constant negative Gaussian curvature) Fig.13b, and show a mixed effect
on the torus that has an external and inner regions of respectively positive
and negative Gaussian curvature, Fig.13c. On spheres, branching patterns are
more bent towards each other as the sphere radius decreases (and the curvature
augments). On the pseudo-sphere, the same tree, represented at three different
altitudes, shows very contrasted shapes. As the pseudo-sphere is of constant
curvature, these variations are not due to a change in the local intrinsic
geometry (that keeps constant everywhere), but due to the change in the
extrinsic geometric component only (indeed, while their product is constant,
the principal radii of curvature are not constant throughout the pseudo-sphere
surface).

Figure 13: Growing tree structures on curved surfaces. (a-c) trees created
with a shooting algorithm to solve an IVP (F primitive). (a) Tree on spheres
(constant Gaussian positive curvature) with decreasing radii (1, 0.7, 0.5).
(b) Tree on a pseudo-sphere surface (constant Gaussian negative curvature) at
different altitudes showing the effect of a local change of the extrinsic
geometry on tree structures. (c) Tree growth on a torus. Left: reference tree
grown in a flat space. Middle: the tree trunk is aligned along the external
great circle (region of positive Gaussian curvature). Remark in the central
region at the tip that the small branches form a very densely organized fan.
Compare with Right: the tree trunk is aligned along the inner great circle
(region of negative Gaussian curvature). In the central region at the tip that
the small branches form a less dense fan. (d) Tree representing the veins of a
leaf, created by joining pre-specified (red) points on the rim (left) to a
main branching system using a RiemannianLineTo primitive (BVP). Next to right:
resulting branching system in the same view as left, followed on the right by
a view slightly tilted, and to the right-end, the back of the leaf.

Branching patterns can also be produced by joining the current position of the
turtle to target points during the construction of a central stem, Fig.13d.
For this the turtle must be placed at an initial position on the surface,
pointing in an initial direction (corresponding roughly to the direction of
the future main stem). Assume that a set of $N$ target points is defined by
their $(u_{n},v_{n})$ coordinates on the surface. We aim at creating a main
stem composed of $N$ segments that follow a geodesic of the surface and such
that, at the end of each segment $n$, a lateral segment is drawn to the target
point $(u_{n},v_{n})$, (List.10):

⬇

1N = 10 # Number of target points

2target_pts = [[0.,0.3],[0.4,0.],...] # Array of (u,v) coords of target points

3ilen = 0.1 # Length of a segment between two branches

4

5Axiom:

6 nproduce SetSpace(Patch(leafblade)) # Nurbs patch in the form of a leaf
blade

7 nproduce InitTurtle([0.,0.,1.,0.]) # initial (u,v,p,q) coords of the turtle

8 nproduce A(0)

9

10derivation length: N+1 # Number of segments on main stem

11production:

12A(n) :

13 if n<N:

14 nproduce F(ilen)

15 nproduce [RiemannLineTo(target_pts[n],20)] # segment made od 20 subsegments

16 nproduce A(n+1)

Listing 10: Tree patterns based on BVP problem (see Fig.13d)

Segments on the main stem are iteratively computed as a series of IVP using a
F statement (line 14). The length of the each segment is determined by the
user. At each step $n$, a lateral segment is computed as a BVP to the $n$-th
target point using the RiemannianLineTo primitive (line 15). This determines
automatically all the points of the lateral segment, and in particular a
specific insertion angle on the main stem (log map, see above), Fig.13d.

###  4.3 Applications

Let us now illustrate on some examples how Riemannian L-systems can be applied
to the modeling of some biological organisms and patterns.

##### Filamentous growth.

Pollen grains are transported from flowers to flowers by wind or animals. They
can germinate if they land on specific elongated cells, called papillae
located at the tip of the stigma, the female organ of the flower. When they
germinate, a pollen tube starts to grow out downward the papillae, and keeping
at the papillae surface [64]. Papillae have roughly a pin-like structure, but
may vary in shape within or between species and present either convex or non-
convex forms. Biologists try to understand the possible physical or chemical
clues that guide the growth of the pollen-tube downwards. One of the
hypothesis is that the precise geometry of the papillae may play an important
role in the guidance of the tube and that the tube could follow geodesics of
the papillae surface [63].

Riemannian L-systems can be used to analyze this growth process and for
example to explore how the shape of the papillae impacts the possible geodesic
trajectories that a tube might follow. Let us consider for example a pin-like
structure and decrease progressively its neck (i.e. its diameter at mid-
height) to pass from a convex to a marked non-convex shape, Fig.14. Starting
from a point $P$ at the tip of the structure representing the position of the
pollen grain and a given orientation of the initial germination of the tube
from the grain, we observe that the trajectory of geodesics is highly affected
by the change in the surface shape: as the neck narrows down, the geodesic
increasingly spirals when approaching the neck and continuously goes downward.
Surprisingly, below a certain neck diameter, the geodesic is reflected and
coils back without crossing the neck any longer, Fig.14 last two right
examples. This is again a consequence of the Clairaut theorem, Equ.35.
Altogether, this suggests that geometry is potentially a key, genetically
regulated, actor in the driving of the pollen tube. This hypothesis was
explored in a recent work using these tools, see [63]

Figure 14: Geodesics on a pin-formed surface. A geodesic of constant length
$l$ in yellow is initiated at a constant angle (28 degrees downward) with
respect to the latitude at the level of the pollen grain position (green
point). From left to right: the mid-height neck of the pin surface is
progressively reduced from left to right (the value of $D$ decreases). As a
result the geodesic coils increasingly around the neck, up to a point where it
cannot pass the neck anymore and coils in the top region.

##### Branching system growth.

Leaves vascular networks are essential to gaz, water and sugar exchanges in
plants. They are an integral part of the development of the leaf and result
from developmental mechanisms that are not yet well understood. In the last
decade however, a few studies have made progress on the understanding of the
genetic regulation of leaf development and its connection with the
construction of the vascular network [67, 31]. In these works, leaves are
considered as flat medium represented by a triangulated mesh in which the
vascular network is embedded. From a geometrical point of view, "the modeling
of leaves that are curved as well as lobed or serrated remains an open
problem. […]. Challenges are posed by […] the need to replace straight
segments with their more complicated counterparts defined on curved surfaces:
the geodesic curves" [67]. Riemannian L-system provide a first generic tool to
address this problem with high-level programming constructs. We illustrate
this on the vein network of the curved leaf of white cabbage, Fig.15a.

For this, a geometric model of the curved leaf blade is first constructed
using NURBS parametric surfaces (also called NURBS patches, [53]) embedded in
the software platform L-Py [6], Fig.15b. Then, the turtle is positioned at the
bottom of the leaf, with its head vector $\mathbf{H}$ pointing roughly in the
direction of the ridge of the surface (lines 4-7 in the code below).

⬇

1N = 10 # Number of branches on main stem

2dl = 0.2 # segment length between 2 branches

3Axiom:

4 nproduce SetSpace(leaf_patch) # space that forms the leaf blade

5 nproduce PlotSpace() # Plots the NURBS patch

6 nproduce InitTurtle([0.001,0.27,1,0.45])

7 nproduce -(6.5) A(N) # correct slightly turtle’s head

8

9# Main apex

10A(n) # n=seg cpt-down on this axis

11 if n > 0:

12 r = BASERADIUS*n/N

13 nproduce F(dl,r)

14 d = 8 + 0.1*n**2 # deflection angle from geodesic @next order

15 a = 5*n/N # insertion angle adjustment

16 nproduce [+(30+a)B(5,r,d) ]

17 nproduce [-(30) B(5,r,-d)]

18 nproduce A(n-1)

19

20# Branch apices

21B(n,r,d): # n=seg cpt-down,r=radius, d=deflexion angle

22 if n > 0:

23 nproduce +(d) F(dl,r*n/N)

24 if n == 4 : # the second segment forks

25 a = -8 if d < 0 else 8 # new deflection angle

26 nproduce [+(5)B(n-1,r,a)]

27 nproduce [-(7)B(n-1,r,a)]

28 else:

29 nproduce B(n-1,r,d)

Listing 11: Cabbage leaf as an IVP (see Fig.15c-d)

A main apex produces segments of equal size following a geodesics (i.e. using
a F primitive) and then producing two opposite lateral branches with insertion
angle close to 30 degrees (lines 13-18). Each lateral branch is composed of 5
segments, the second of which forks into two branches (lines 23-29).

The different parameters of this model have been adjusted by hand to visually
reproduce main geometric traits of the image in Fig.15a. However, it can be
noticed that the number of parameters remains limited and they only
necessitate fine grain tuning (at the level of branch insertion angles in
particular). The core branching system makes use of geodesic lines and is
straightforward to program. Only deflections with respect to this geodesic
pattern need to be adjusted. In particular, the main stem was not adjusted at
all (line 13) and the corresponding geodesic follows naturally the ridge of
the curved leaf blade, despite the fact that it is slightly twisted near the
tip of the leaf as can be observed in the top of Fig.15d.

Figure 15: Cabbage leaf model. (a) photos of a white cabbage leaf (up: top
view, below: side view). (b) Approximated NURBS model of the cabbage leaf
(c-d) different views of the vascular network constructed with Riemannian
L-systems, with veins corresponding to geodesics computed as IVPs.

##### Plant branching system.

Some plants, such as lianas or ivy, cannot support their own weight as they
grow. They then often grow on a supporting structure, which could be another
plant, a rock or a building for instance. The supporting structure, in turn,
may possess intricate geometric characteristics. If we assume that, in the
absence of any extra external force the axes of the plants grow straight, then
we can model the growth of the axes by following the geodesics of the
supporting structure. Conventional approaches rely on collision detection
techniques using a voxel representation of the 3D space [23] or bounding
volume hierarchies [74] to determine the path on the structure that should be
followed. In the case of complex geometry of the supporting structure,
substantial refinements of the voxel representation are required to capture
the detailed changes in curvature leading to computationally intensive
generation. Riemannian L-systems present an efficient means to simplify such
modeling process by embedding the generated plant shape directly onto the
supporting structure, thus mitigating the complexities associated with
intricate geometry of the supporting structure.

In the example depicted in Fig.16, the climbing of an ivy plant is directly
simulated on the external surface of a tree trunk. This process involves
representing the trunk as a generalized cylinder characterized by an S-shaped
central axis, and with a linearly decreasing radius along the axis. The
cylinder cross-section incorporates concavities indicative of the differential
radial growth of the trunk and its fusion with aerial components of the roots.
The initiation of the ivy growth occurs at the base of the trunk, with an
upward orientation. The simulated trajectory then follows naturally the ridge
along the trunk surface. Ramifications are regularly generated at a constant
insertion angle along the axes, resulting in the coverage of the trunk surface
by the ivy structure.

Interestingly, the disposition of ivy leaves is not confined to the surface of
the shape. For this, the embedding of the generated shape onto the surface can
be stopped using the EndSpace command. When incorporated into the L-string,
the geometric interpretation of the subsequent modules of the string will
recover the conventional behavior of L-systems within the euclidean 3D space
embedding. For the ivy example, the normal vector to the surface at the leaf
insertion point is estimated, and the turtle is reoriented accordingly.
Subsequently, a roll rotation is then executed to establish a horizontal
orientation for the leaves. The petioles and the leaf blades are then
generated. This example illustrates how multiple embeddings can be combined
simply within the same Riemannian L-system. This could be generalized further
by considering the generation of a plant onto multiple support structures. For
this, the different embedding spaces can be set by inserting several SetSpace
or EndSpace commands into the L-string.

Figure 16: Climbing ivy. A tree trunk is modeled as a generalized cylinder on
which the growth of an ivy is simulated. On the left, the wireframe
representation of the trunk and the branching system of the ivy. On the right,
semi-transparent polygonal representation of the trunk with the full leafy ivy
structure.

##  5 Feedbacks between surface and embedded forms

Be they of mathematical, physical or biological nature, forms are built
according to construction rules that operate within some substrate space, e.g.
a plane or a curved surface. These rules may be completely independent of the
space in which the construction process operates. This was the case in all the
previous examples. However, in many situations, morphogenesis relies on
external cues that are used during development by the construction rules to
make decisions and orient form development (by changing growth orientation or
speed, by creating branches, etc.). In particular, in such cases, geodesics
can be seen as default trajectories and construction rules indicate how to
deviate from these reference trajectories to construct the target form based
on local information.

Spatial informational cues may be considered as fields living in the substrate
space, and that may be constant or change in time. In general, these fields
have scalar, vectorial or more generally tensorial types, and may represent
either geometric (e.g. curvature, principal directions), physical (e.g,
obstacles, molecular concentrations, stresses, material properties) or more
abstract (e.g. energy density, directional anisotropy, signals) local
quantities.

Forms thus, in general, result from the interaction between three key factors:
A substrate space on or in which the form develops, information fields defined
on this space that may or not be (partially) produced by the form itself, and
the growing form (here a filament or a branching network). Hereafter we
explore how sensing external fields may contribute to shaping forms in
Riemannian L-systems, first in a fixed embedding space, then in a dynamically
deforming embedding space.

###  5.1 Information feedback

In Riemannian L-systems, to make use of external fields living in the
substrate space (here a surface) during development, we must allow the form
being constructed to probe its space environment at anytime. This will be
carried out by allowing the modeler to access turtle geometric information
from within the production rules.

For this we extended the L-system mechanism of sensitivity to the environment
(see section 2.3). By using a new special query module ?T inserted in the
L-string constructed at a given derivation step, the modeler can ask the
turtle to store its complete current state during the interpretation phase
(like for the modules ?P, ?H, ?U used in context sensitive systems - see
section 2.3) in this query module. At the next derivation step, the modeler
can then access the turtle’s state recorded in the variable of module ?T and
take decisions based on the information contained in this state.

At any moment, the turtle variable stored in ?T contains information about the
turtle state itself and about the current space (i.e. the surface) in which it
moves. In this way, the space can be queried to retrieve all types of
information attached to the surface: geometric primitives (covariant basis,
surface normal, principal directions, Gaussian and mean curvature, principal
curvature and directions …) using dedicated primitives, as well as any type of
field value stored by external processes on the surface.

Figure 17: Making use of information available in the embedding space. (a-c)
Random walks in regions of positive Gaussian curvature. (a) Curved space with
seven random walkers initially positioned at the indicated frames. (b)
Canalyzed random walks after 200 steps for each walker. (c) Map of Gaussian
curvatures $K$: red regions with $K>\epsilon$, and white to blue for
$K\leq\epsilon$, $\epsilon=0.002$. (d-f) Self-avoiding branching structures
(d) on a sphere, (e) filling the sphere, (f) on an egg-box like landscape.
Colored points represent growing apices

##### Scalar fields

Scalar fields can express either local geometric properties of a surface or
spatial distributions of physical or biological quantities defined on the
surface. Let us consider for example how the random movement of a set of
turtles can be canalized by their ability to read out locally the geometric
characteristics of the embedding surface and to use it to make move decisions,
Fig. 17a-b. We assume that the turtles locally sense the surface curvature and
avoid to go in flat (0 Gaussian curvature) or saddle-like regions (negative
Gaussian curvature). These random walks are thus restricted to regions of
(strictly) positive curvature (indicated in red in 17c). To model this
situation in Riemannian L-systems, one would basically modify the random walk
algorithm described in List.2 as follows:

⬇

1epsilon = 0.002

2production:

3?T(state) A(n) :

4 u,v,p,q = state.uvpq # Current turtle’s position (u,v) & orientation (p,q)

5 surface = state.space

6 found = False

7 while not found

8 a = 360*random()

9 uvpq_rot = riemannian_turtle_turn(t.uvpq, surface, a)

10 uvpq_seq = riemannian_turtle_move_forward(t.uvpq, surface,step)

11 uvpq_new = uvpq_seq[-1] # Last point of the tentative move sequence

12 K,H,kmin,kmax = surface.localCurvatures(uvpq_new[0],uvpq_new[1])

13 if K > epsilon: found = True # Test if Gaussian curvature is positive

14 nproduce +(a)F(uvpqs) ?T A(n+1)

Listing 12: random walk keeping on regions of positive Gaussian curvature (see
Fig.17a-c)

The query module ?T is inserted just before each moving apex in the L-string
at each derivation (line 14). At the next derivation, the turtle’s state (line
3) automatically updated at the last interpretation step can be used to get
the current position and orientation of the turtle (line 4) as well as the
current embedding surface (line 5). Then the model tries to make moves in
different random position until one of this moves happens to be in a region of
positive curvature (line 7-13). Different primitives make it possible to
instruct the turtle to turn by an angle $a$ (line 9), to return the sequence
of coordinate corresponding to a potential step in this new direction (line
10), to compute the the Gaussian curvature $K$ at the destination point (lines
11-12) and to test whether its value is positive (line 13). As a result, one
can observe that the random walks get canalized in the regions of positive
Gaussian curvature, without being able to cross areas of flat or negative
curvature Fig.17b-c.

A scalar information can also result from a read out of the physical
environment. A geometric form for example can be physically constrained in its
development by its already existing parts. As an example, consider the
construction of self-avoiding branching patterns, Fig.17d-f. A specialized
data-structure and associated primitives have been developed to record
segments produced as the form grows and to test the presence/absence of
already existing segments on the surface at specific locations, List.13. An
empty data-structure is first created (line 10). Then the apices of branching
structure grow (line 16). Like before, a query module is used to recover the
current state of the turtle before each apex. Instead of creating a segment,
an apex first computes the path on the surface that would correspond to this
new segment (line 17). It then test whether an intersection is detected with
previously created segments (line 18). If yes the apex does not grow (line
19). If no, the new segment is added to the list of already existing segments
(line 21), and the segment is drawn (line 22), together with a new lateral bud
(line 24) and a renewed main apex (line 25).

⬇

1R = 1. # Radius od the Sphere

2N = 7 # Depth of the tree recursion

3iangle = 45 # Insertion angle

4ilen = 0.2 # Length of a segment between two branches

5trajectories = None # Will contain the set of already existing tree branches

6

7Axiom:

8 space = Sphere(R)

9 nproduce SetSpace(space)

10 trajectories = LineSet(space) # Initializes the set of existing branches
(lines)

11 nproduce InitTurtle([0.,0.,0.,1.])

12 nproduce A(0)

13

14derivation length: N

15production:

16?T(ts) A(n) :

17 uvpq_s = forward(tstate,ilen) # Pre-computes a new segment (sequence of
points)

18 if trajectories.test_intersection(uvpq_s): # Does the new segment intersect
previous ones?

19 nproduce ?T A(n+1)

20 else:

21 line_id = trajectories.add_line_from_point(ts.uvpq,uvpq_s) # updates set of
trajectories

22 nproduce F(uvpq_s) # Draw the precomputed new segment sequence of points

23 a = iangle if n % 2 else -iangle

24 nproduce [+(a) ?T A(n+1)]

25 nproduce ?T A(n+1)

Listing 13: Non self-intersecting tree (see Fig.17d-f)

This procedure is generic, and does not depend on the underlying space.
Fig.17f for example illustrates the same self-avoiding branching algorithm
applied to an egg-box surface with a more complex geometry with both positive
and negative Gaussian curvatures.

##### Vector fields

Scalar fields provide information that can indirectly be used to affect
turtle’s trajectory. A more explicit directional information may be provided
by a vector field. At any point of the surface the turtle can compute the
angle $\alpha$ between its current heading direction $\mathbf{H}$ and the
local vector value of the field $\mathbf{V}$, and orient its trajectory
locally in the direction indicated by the vector field. This deflection of the
otherwise geodesic trajectory in the direction of a vector field is in general
called a tropism. Formally, let us denote $\alpha$ the current angle between
the turtle heading direction $\mathbf{H}$ and the local value of the vector
field $\mathbf{V}$ at the turtle’s position,
$\alpha=\widehat{\mathbf{H},\mathbf{V}}$, and $\Delta s$ the small step length
that the turtle must make at the next derivation step. We assume that the
tropism will deflect the turtle’s initial direction from a geodesic by
imposing a geodesic curvature proportional to the angle $\alpha$ and that
$\Delta s$ is sufficiently small so that the geodesic curvature is considered
constant over the step length $\Delta s$. Then, we have:

$\frac{\Delta\alpha}{\Delta s}=\sigma\alpha,$ (36)

where $\sigma$ is a sensitivity to the "tropism force". This provides a direct
expression of the deflection angle $\Delta\alpha$ by which the turtle’s
movement must be affected at the next iteration step to progress over a
distance $\Delta s$ due to tropism.

To illustrate this, consider a geodesic starting at the equator of an
ellipsoid of revolution and heading east and slightly north, Fig.18a. In the
absence of a field (or of a feedback between a field and the turtle’s
displacement), the turtle follows a geodesic trajectory Fig.18a. In the
presence of a field (Fig.18b), the turtle may react to the field and use this
information to locally modify its natural trajectory. Here, we consider a
vector field resulting from the gradient of a scalar field corresponding for
example to the diffusion of a substance from the north pole (red, resp.
yellow) is high (resp. low) concentration of the substance). By reading this
gradient the turtle can deflect its trajectory in the corresponding direction.
As a result the trajectory is attracted at the north pole and trapped in a
circular attractor (Fig.18b). Interestingly, this attractor results from an
equilibrium between the effect of the surface geometry that, in the absence of
tropism ’force’, would keep the turtle on a geodesic and thus make it possible
to escape the north region after having visited it (Fig.18a), and the tropism
force that constantly deflects the trajectory towards the north pole.

Figure 18: Deflection of geodesic trajectories using surface vector fields.
(a-b) Tropism on an ellipsoid of revolution (a) No field: the trajectory
(green) is a geodesic starting at the equator and heading east, bending 30
degrees north. (b) presence of a scalar field (red = high, yellow = low
values). The trajectory, with identical initial conditions converges to a
circle in the north region. (c-d) Simulated tip-growing filament trajectories
on pin-formed structure with a scalar field at the surface (color gradient
from red (high) to dark blue (low values)). (c) left-most: Geodesic trajectory
(not interfering with the scalar field). To the right: effect of tropism
resulting of an interaction with the scalar field. Geodesics are deflected in
the direction of the gradient of the scalar field, with an increasing
intensity $\sigma$ from 0 to 5. (d) The wireframe structures in the bottom row
show whole trajectories.

Programming such a feedback function follows the design pattern illustrated in
List.14. A query module ?T before the growing apex makes it possible to
retrieve at each step the current state of the turtle at the apex (lines 16
and 27). At the position $(u,v)$ of the apex, the value of the tropism vector
field is evaluated (here as the gradient of a scalar field, lines 18-20). The
angle $\alpha$ between the tropism vector and the heading direction of the
turtle is then computed (lines 22-23) as well as the deflection angle due to
tropism over a step length $slen$ according to the tropism attraction
described in Equ. 36 (line 24). The sign of this angle depends on the relative
orientations of the heading and tropism vectors (line 25). Finally, the turtle
turns according to the deflection angle and draws a portion of geodesic over
length $slen$ (line 27), thus deflecting its trajectory towards the tropism
vector.

⬇

1a, b = 1., 0.5

2N = 200

3sigma = 1. # sensitivity to gradient

4angle0 = 45 # Initial angle

5slen = 0.1 # Length of a step at each derivation

6def field(u,v): # function returns a scalar value as function of u,v

7 ... # \--> not detailed

8

9Axiom:

10 nproduce SetSpace(Ellipsoid(a,b))

11 nproduce InitTurtle([0.,0.,1.,0.])

12 nproduce ?T A(field)

13

14derivation length: N

15production:

16?T(state) A(sf) :

17 u,v,p,q = state.uvpq # retreives current pos, dir of turtle

18 gradf_uv = gradient(state,u,v) # gradient of the scalar field at u,v

19 pushfwd = ts.pushforward(u,v) # local pushforward operator

20 gradf = pushfwd(gradf_uv) # gradient vector on the surface in 3D

21

22 heading = state.heading # turtle’s head in 3D

23 alpha = angle_between(gradf, heading) # deviation of turtle’s head from
gradient

24 defl_angle = sigma * alpha * slen # deflection angle to correct turtle’s
dir

25 sgn = 1 if state.space.positive_orientation(heading,gradf)) else -1

26

27 nproduce +(sgn*to_deg(defl_angle) F(slen) ?T A(sf)

Listing 14: Tropism on ellipsoid of revolution (see Fig.18a-b)

This design pattern can be applied to model similar situations in biological
applications. For example, (Fig.18c-d) illustrates the possible effect of a
surface tropism on a tip-growing filament at the surface of a pin-like
surface. This could for instance represent the growth of a pollen tube that
would be attracted by a chemical gradient. In the absence of tropism, the
filament follows the geodesics of the pin-like surface Fig.18c,d-left. In the
presence of a tropism due to the surface gradient (indicated by a variation of
colors from red to blue), the filament is deflected from geodesic curves and
heads more rapidly towards the bottom of the surface. This phenomenon is
amplified by an increase of the coupling between the filament growth and the
directional cue (increase of $\sigma$ from left to right in Fig.18c,d).

###  5.2 Growing forms on dynamic surfaces

Forms can be constructed on surfaces that change in time. This change may
affect the embedded form in different ways that we explore in this section.
Here we restrict our analysis to form construction processes that are
significantly more rapid than the change of the surface geometry so that,
during the form construction, the surface geometry can be considered steady.
In this way, after one growth step of the surface, several steps of growth for
the form can be simulated repeatedly until convergence.

##### Convected versus floating forms

Once a form has been constructed on a surface, two extreme subsequent types of
evolution of this form can be considered.

Let us first assume that the points on the surface can be tracked in time.
They can be considered as material points whose movement define the surface
evolution in the 3D space. At the initial time, the form passes through a set
of material points to which the form is considered to be "attached". While the
surface changes in time, the relative positions of the material points change,
which induces a corresponding deformation of the form geometry on the surface.
As the form depends on the definition of material points and their movements,
we say that such form definition is extrinsic. An extrinsic form is being
convected by the material flow on the deforming surface.

This situation is illustrated on Fig.19a-b using a fractal form. A Sierpinski
carpet, e.g. [52, p 81], whose construction process is depicted in Fig.19a, is
embedded in an initially planar surface that will progressively take the form
of an egg-box. This can be modeled by allowing the parameters of a parametric
surface, here a NURBS patch $x^{i}(u,v,\lambda(t))$, where $\lambda(t)$
represents the surface parameters at time $t$, to change smoothly in time, and
to consider that the uv-coordinates correspond to coordinates of material
points. We can observe that the deformation of the Sierpinski carpet does not
change its topological structure, and only its geometry is smoothly affected,
reflecting the smooth deformation of the embedding surface Fig.19b. In
Riemannian L-systems, this form convection can be achieved using indirect
interpretation of L-strings (see section 3.4). In this mode, the form is not
directly drawn as a sequence of geodesics on the surface. Rather, it is drawn
in the uv-parameter space and then projected on the surface, reminiscent of
texture mapping in computer graphics, e.g. [14]. In this way, the form drawn
in the uv-coordinates is that of Fig.19a. The line segments projected on the
surface are in general not geodesics, but the projection preserves material
points neighborhood. As a result, the form on the surface appears simply to be
convected by the smooth surface material deformation that preserves its
integrity and topology. To keep track of the uv-coordinates of a previous
drawing, we introduced a primitive StaticF that behaves like a classical
forward instruction F(x), i.e. it computes a geodesic of length $x$ from the
current position and heading of the turtle. However, the first time it is
called, it memorizes the uv-positions of the computed geodesic. Then, at next
derivation steps, on further calls, it will not recompute a geodesic on the
modified surface. Instead it will use the previously cached geodesic
information.

Figure 19: Feedback of surface dynamics on forms constructed at the surface.
(a) Steps of the construction of a Sierpinski carpet in the Euclidean plane
corresponding constructed with a L-system. At each step, the pre-fractal form
is obtained by the trajectory of a turtle moving from A to B along convoluted
paths. At scale 0 (left), the form is approximated by a simple segment. The
turtle draws this segment by going straight from point A to point B. Then, at
scale 1, this segment is refined into 8 smaller segments of length 1/3 of the
original segment length each, as illustrated on the next diagram. The turtle
draws this pattern by following the trajectory indicated by the grey arrows.
The refinement process then continues at higher scales by decomposing further