[Nature of Code] Part 2: Fractals and Cellular Automata
Nature of Code note series
Fractals
The term fractal is coined by mathematician Benoit Mandelbrot in 1975. It means broken or fractured.
A structure needs to have two properties to be a fractal:
- Self-similarity: when you zoom in on it, it looks similar to (or exactly the same as) the whole.
- Fine structure at small scale.
- Generated from a recursive definition.
A binary tree with infinite depth is an example of exact fractal.
The stock price curve is an example of stochastic fractal. When you zoom in, the shape isn’t exactly the same, but it has the same quality, a quality of squiggly randomness. The coast line is another example.
A famous and extremely simple fractal pattern is the Cantor Set. It divides a line into 3 equal segments and remove the middle one. You get the Cantor set by doing this to every segment recursively.
Koch curve: the “monster curve”
The generation rule of the Koch curve is very similar to the Cantor set:
- Divide a line into 3 equal segments
- Remove the middle segment
- Replace the middle segment with an equilateral triangle
- Repeat this to all the line segments.
In code, the way to do this is to use two lists, current
and next
. Pseudocode:
Start with current = [line0], next = []
i = 0
while i < max_steps:
for line in current:
find all the points of interest on line: a, b, c, d, e
(a is the left end, b is at the 1/3 mark, c is the tip of the triangle,
d is the 2/3 mark, e is the right end)
create the new line segments (ab, bc, cd, de) and add them to next
current = next
next = []
i++
Binary fractal tree
Generation rule:
- Start with a vertical line, go to its tip
- At the tip in 1, rotate 45 degrees and draw a new line with half the length
- At the tip in 1, rotate -45 degrees and draw a new line with half the length
- Go to the tips of 2 and 3, repeat.
In pseudocode,
branch(next_length):
drawLine(0, 0, 0, -next_length)
translate(0, -next_length)
len = next_length/2
if len > threshold:
pushMatrix()
rotate(theta)
branch(len)
popMatrix()
pushMatrix()
rotate(-theta)
branch(len)
popMatrix()
In each recursive step we call branch()
.
pushMatrix()
and popMatrix()
Review: Notice that pushMatrix()
creates a new reference frame that saves the current transformation (the current translate and rotate state).
To visualize a reference frame state, it can be represented by a unit vector which has the starting point at the new (0, 0)
and pointing to the new (0, 1)
direction.
Then popMatrix()
can forget about all the transformations performed after the last pushMatrix()
, and restore to the previous state. Thus it’s a stack of transformation “sandboxes”.
These are Processing terminology. In P5.js they are called push()
and pop()
.
L-Systems: applying recursion to text
An L-system or Lindenmayer system is a parallel rewriting system and a type of formal grammar. Lindenmayer introduced L-systems in 1968 to describe the behaviour of plant cells and to model the growth processes of plant development. L-systems have also been used to model the morphology of a variety of organisms and can be used to generate self-similar fractals.
An L-system has 3 components:
- Alphabet: characters allowed in this L-system
- Axiom: an initial string
- Rule set: a mapping of character to its descendant(s)
For instance, we have an L-system that has
- Alphabet: A, B
- Axiom: A
- Rule set: A -> ABA, B -> BBB
How does it work?
Gen 0: A
Gen 1: ABA
Gen 2: ABABBBABA
Gen 3: ABABBBABABBBBBBBBBABABBBABA
...
If we think about it as drawing things, set A to be drawing a line -
, and B to be moving forward _
(whitespace), then what is the effect above? Let’s see,
Gen 0: -
Gen 1: - -
Gen 2: - - - -
Gen 3: - - - - - - - -
...
We have recreated the Cantor set!!
This means that we can design L-systems that create fractal patterns.
Space Colonization
Space Colonization is an algorithm that creates a fractal tree that grows by capturing pre-generated “leaves”. These leaves are usually generated randomly to fill up the space, and they attract the branches to grow towards them.
Here’s the pseudocode of this algorithm for 2D space.
class Leaf:
constructor():
Create leaf with random 2D coordinates (x, y), and a boolean flag
`reached`.
show(): Render the leaf on canvas
class Branch:
constructor(parent, position, direction):
A branch is represented by its position, its parent branch's position,
and a direction copied from its parent. This direction can be used to
create new directions for its children.
Properties of a Branch object:
parent, pos, dir, originDir (keep track of parent's direction),
count (# leaves that find this branch as closest),
stepSize (a parameter to tune to grow faster or slower)
reset(): Set the direction back to originDir (parent's direction),
and count to 0.
show(): If parent exists, render the line from parent position
to position
class Tree:
constructor():
A Tree has properties: an array of leaves, # leaves, an array of
branches, a maxDistance, and a minDistance.
It needs to be initialized. The initialization contains several steps:
- It needs a root as the 1st element in branches. This root is a Branch
object with a position, direction and no parent.
- All leaves need to be generated in a loop.
- As a prerequisite of grow(), first the root needs to grow to the
vicinity of some leaves in order for the leaves to attract them later.
This process can use a while loop and a boolean variable `found`
initialized to false. Set current branch to root. While not found:
Loop through all leaves and check the distance to
the current branch, if it's less than maxDistance, set found to
true. If after all leaves, found is still false, grow a new branch
on the current branch that has the same direction as current branch,
and set this new branch to be the current branch.
This way we can grow the root to the vicinity of some leaves and keeping
its initial direction unchanged.
grow():
This is the key part of this space colonization algorithm.
Loop through all leaves. For each leaf, loop through all branches and
find its closest branch that is within the (minDistance, maxDistance)
interval. This closestBranch is the one this particular leaf attracts.
In the double for-loop, calculate the distance between the current pair
of branch and leaf. If it's less than minDistance, set leaf.reached to
true, and closestBranch to null. Else, compare this distance to the
previous record min distance, if (it's lower or closestBranch is null),
set closestBranch to current branch and update the record.
After looping over all branches for a leaf, IT IS THE PART THAT
GENERATES THE NEW BRANCH:
- Get the vector (closestBranch -> leaf), normalize it
- Add this normalized vector to closestBranch's direction which will be
used as the direction of its child branch. Note that this direction
is cumulative with all previous or future leaves that share this
closestBranch.
- Increment clostBranch's count which keeps track of the number of
leaves that share this closestBranch, in other words, influenced the
cumulative direction.
After the double for-loop that runs through all branches for all leaves,
check all leaves with one more pass. If it's reached, remove it from
the array. (JS Trick: to remove element while iterating on the same array,
iterate backward)
Next, iterate over all branches. If a branch has count > 0, make a
child branch from it using its updated direction. One trick to make the
tree better looking is to shorten the length of a new branch if it's
shared by a lot of leaves, by dividing its length by its count. Then,
don't forget to reset the branch!
show(): Call show() of all leaves and branches.
To help this tree reach more leaves, we can add a small random perturbation in direction whenever we grow a new branch.
You can check how this works in action here on my website.
A Variation of Space Colonization
I find it satisfying to look at perpendicular lines in 3D space so I tweaked the above algorithm a bit to grow a fractal tree with branches in orthogonal directions, x, y and z. The place I tweaked is the new branch growing part.
The idea is simple. Here’s the pseudocode for the new branch generation part.
- Get the vector (closestBranch -> leaf) as before, normalize it.
- Find the x, y, z component that has the largest magnitude, and set it to be
the new direction
That’s it. The way to make it 3D is even simpler. Just add a z
component to all geometry.
Check out the result here. I added the logic that once all branches are grown, skip the grow()
method in the draw loop. However, since there are many branches in 3D space, the performance is still not very good. If you have suggestions to improve this, please let me know.
Cellular Automata
The most famous cellular automata algorithm is the Game of Life.
A cellular automaton has a grid of cells, each cell has
- a state: dead (0) or alive (1)
- a neighborhood: adjacent cells
The state of a cell at time t
depends on its neighbors in t-1
. So we have the concept of generations.
Wolfram classifies the possible outcomes of a cellular automaton into 4 types:
- uniformity
- oscillation (repetition)
- random
- complexity
With certains rules (such as rule 30), even a 1D Game of Life can be a pseudorandom number generator.
You can check out the artilce that lists the rules on Wolfram’s website. Some rules are particularly interesting. For example, rule 90 generates a Sierpinski Triangle which is a fractal!
Game of Life (2D)
There are 3 scenarios in Game of Life
- Death
- Overpopulation
- Loneliness
- Birth
- Exactly 3 live neighbors
- Stasis
- In all other cases, stay the same
For the edges, we can either ignore them or wrap them around, i.e. treat position -1
as position length
for each row.
If you think about it, every digital image processing algorithm is a cellular automaton! Each pixel is essentially a cell. A convolution, for example, is a cellular automaton! So is a Gaussian blur, a ripple effect, and so on.
An interesting expansion to the simple black and white Game of Life is to make each cell an object and let it store its history. We can make it blue if it’s within N steps of birth, red if it’s within N steps of death, for instance.
Some more possibilities
- The cell could be non-rectangular, could be hexagonal, triangular, etc.
- Add probability, e.g. a rule of 80% of dying.
- Continuous states
- Image processing: an application to Cellular Automata, say, a water ripple effect.
- State history
- Moving cells. Change movement based on neighboring cells (flocking)
- Nested complex systems, say a cell is another CA.
Reference
- The Coding Train videos