Generating a 3D growing tree using a space colonization algorithm

Steps of generation

In this article we will create a procedural growing tree using a space colonization algorithm. The article will be divided into 2 parts: the space colonization algorithm and the generation of a 3D mesh based on the data provided by the colonization. The implementation will be done in C# within the Unity Engine. Source could is available at the end of the article.


1. Generating a tree: the space colonization algorithm

There are many algorithms to generate realistic trees. Daniel Shiffman covers two of them (fractal, L-system) in The Nature of Code, Chapter 8 [2]. In this article, we will use a space colonization algorithm. The algorithm I will walk you through comes from a paper, Modeling Trees with a Space Colonization Algorithm [1]. If reading a scientific paper isn’t a burden for you, I suggest that you read it first and then come back to read my approach, even though this step isn’t mandatory.

Figure 1. The different steps of the space colonization algorithm

The space colonization algorithm can be described as following:

  1. A volume inside which the leaves will be contained is defined. Its shape will have a direct influence on the shape of the final tree
  2. This volume is populated by leaves (points in the space), using any kind of distribution
  3. A branch starts to grow from the bottom of the volume towards the leafs
  4. The branch is attracted by the leaves and starts to split over and over to reach more of them. Once a leaf is reached, it is removed from the volume
  5. This keeps going until no more leaves are left

a. Distribution of the leafs within a volume

For the sake of simplicity, a sphere will be used as a container for our leaves. This makes the distribution easier since a random position within a spherical volume can easily be computed. Other volumes could be used as well as other distributions. Our volume will be defined by 2 constants: its position and its radius.

Figure 2. The spherical containing volume

Leaves are now going to be added within this volume. If you are working with a custom volume, you will have to work your distribution out to fit it properly. In our case, because the volume is a sphere, we can randomize a direction vector and its magnitude. The direction can be computed by generating 2 random angles. Let  \alpha \in [0; \pi] and  \theta \in [0; 2\pi]. Let  \vec{D} be our direction vector:

    \[\vec{D} = \begin{bmatrix}    cos(\theta) \cdot sin(\alpha) \\    sin(\theta) \cdot sin(\alpha) \\    cos(\alpha) \end{bmatrix}\]

Let  r be the radius of our containing sphere. The magnitude  d \in [0; r] is a random value between 0 and  r. Let  P be the random position of our leaf:

    \[P = \vec{D} \cdot d\]

Figure 3. Random distribution of 400 leaves within the sphere
Please note that the mesh of the sphere volume is being drawn after the leaves without any depth test

b. The branches

Before diving into the logic behind the growth of the branches, we need to think about a structure that will suit our needs. Branches can be thought as a segment between 2 points. Moreover, we can already plan what other characteristics a branch should have. First, we should be able to navigate the tree starting by any branch. Thus each branch must store its parent and its children. Because we might need to use its direction quite often, let’s also store it instead of computing it each time we need it. When growing, the branches are attracted by some leaves. We will track those leaves in a list. This is a C# implementation of such a class:

class Branch {
	public Vector3 _start;
	public Vector3 _end;
	public Vector3 _direction;
	public Branch _parent;
	public List‹leBranch› _children = new List‹Branch›();
	public List‹Vector3› _attractors = new List‹Vector3›();

	public Branch (Vector3 start, Vector3 end, Vector3 direction, Branch parent = null) {
		_start = start;
		_end = end;
		_direction = direction;
		_parent = parent;
	}
}

At each iteration, we are going to add branches to the extremity of the existing branches in attraction range of 1 leaf at least. All the branches of our tree will be stored in a list. A first branch is going to be added under the leaves volume (the sphere) and will be used as a starting branch for our colonization algorithm. 2 constants are also going to be set before the initialization: the starting position _startPosition and the length of a branch _branchLength. In this case, we set the position of the volume to [0; 7; 0] and its radius to 5, so the starting position will be [0; 0; 0] and we will define the branch length to 0.2.

In all the following illustrations, the branches will be drawn in green and their extremities in purple.

Figure 4. Our first branch

c. Let our branches grow

This is the core of the algorithm. I will break down this part into smaller steps for a better understanding of the process.

Let’s first forget about the leaves: we want our tree to grow like a bamboo, straight to heaven. At each iteration, we will create a new branch from the extremity, in the same direction as its parent. Since the direction is a vector pointing to the top, our new branch will also point to the top, and its starting point will be the same as its parent end point.

Figure 5. Our first branch growing in the simplest manner

We now want the branches to be attracted by the leaves. To do so, we first need to define an attraction range  A_r. When the distance between a branch and a leaf is lower than the attraction range, the branch will be attracted by the leaf. On each iteration, for each branch, the distance  d_{bl} between the end point of the branch and each leaf is computed. If the extremity is in the attraction range, we will add the leaf as an attraction point to the branch. The following example demonstrates this process only on the last branch added to our tree:

Figure 6. The leaves (in red) are tested to see if they are in the range (in grey) of our only extremity. If so, they are considered as active (in yellow)

For our tree to split in multiple directions, we need to find the attraction points for each branch of the tree. If at least 1 attraction point can be found, a child branch will be created towards the average direction to those points. Let  \overrightarrow{D_{branch_b}} be the direction of the child of  branch_a,  E_{branch_a} the position of the extremity of  branch_a,  P_{attr} the list of attraction points of  branch_a, of size  i_p. The direction of the child of  branch_a can be computed as the normalized sum of the normalized directions between the attraction points and the end of  branch_a:

    \[     \overrightarrow{D_{branch_b}} = \frac{1}{i_p} \cdot \sum_{n=0}^{i_p-1} \frac{ \overrightarrow{P_{attr}[n]} - \overrightarrow{E_{branch_a}} }{ \| \overrightarrow{P_{attr}[n]} -  \overrightarrow{E_{branch_a}}\| } }\]

The next example shows how this operation is applied to one branch, the extremity of our tree. Later on, we will apply this operation to all the branches of the tree, resulting in branches growing in different directions.

Figure 7. The extremity “creates” a branch that grows towards the active leaves

A quick issue arise: when the end point of a branch ends up right on the segment between two leaves in attraction range, the branch can’t grow anymore. Indeed, if you apply the above equation to such a situation, you end up with \overrightarrow{D_{branch_b}} = \vec{o}. To fix this, and to add some overall randomness to the generation, we can add a randomness factor. We will add to \overrightarrow{D_{branch_b}} a normalized random vector multiplied by this randomness factor.

Figure 8. Introduction of randomness into the growth

Okay. Our tree grows indeed towards the leaves, but it wiggles around when it reaches them. Next step needs to be introduced: the kill range. Once a leaf is in the kill range of a branch… well, we kill it. Let’s say, in a more fashioned way, we remove it from the list of attraction points. The kill range  K_r must of course be lower than the attraction range, but should also be greater than the branch length to avoid weird results. This could be illustrated and explained in detail, but such an explanation seems unnecessary.

Figure 8. Leaves are removed when entering the kill range (in dark grey)

We can now extend this process to all the branches instead of only the extremity of the tree. On each iteration, we want each leaf to be attached to the closest branch in its attraction range. Then we create new branches using the same process as before. We also need to cover edge cases: what if no branch is in range of attraction of a leaf ? We can keep tracking of the extremities of the tree. Then if we are in the stated case, we can have those extremities to grow in their current direction with some randomness. This will actually cover the growth of our first branch aswell. We also want to stop the growth is no leaves are left.

Figure 9. First generation of our tree

Oh, that was quite disappointing. Let’s try a new one.

Figure 10. New try

That’s better, but still some issues persist. Let’s point out those issues. First, some leaves can’t be found by our tree. This could be solved by increasing the attraction range and by adding more leaves. We will try to find a great balance for those tweaks. Second, our leaves distribution seems to be wrong. Too much leaves are within the center. The more we are away of the center, the less the leaves density is. This is because the volume of a sphere does not increase linearly along its radius and because we use a linear distribution of the leaves along the radius. We currently pick a random radius  r = Rand(0, 1) * radius. We can find a function between 0 and 1 which returns bigger values the closer it gets to 1. I used a graph tool to find such a function using basic knowledge and by testing:

    \[    f(x) = sin(x \cdot \frac{\pi}{2}) ^ {0.8}\]

Figure 11. Graph of the distribution function

Let’s compare this new distribution to the previous one:

Let’s see how our algorithm performs on 800 leaves with the corrections made to the distribution and attraction range.

Figure 13. Tree generation with an increased number of leaves, better distribution of the leaves and an increased attraction range

d. Optimizations

We can make some optimizations to our algorithm. For instance, we can divide the space into sub-spaces cubes and assign the branches into their corresponding sub-space. When we parse the attraction points, instead of looping through all the branches we can only look for those in the sub-spaces neighbors. I will cover this technique in another article because this can be used in many other cases, and in many different ways.


2. Mesh construction from our branch structure

In this part we will turn our branches graph into a 3D mesh. We will first build the mesh regardless of the size of each branch. Then, we will use our structure to assign a radius to the branches given their distance from the extremities of the tree. And finally, we will animate the extremities so that our tree can grow smoothly over time, while being generated.

a. Static mesh from our branches

First thing: we need our branch data structure to store more informations for our mesh construction. We will store the indices of the vertices in the branch data structure to build the faces later one. The mesh construction is pretty straightforward. For each branch, we will create  S vertices around the end point and along the branch direction. We will then create faces from those vertices. Let  N_b be the number of branches,  N_v the number of vertices of our mesh:

    \[N_v = (N_b+1) \cdot S\]

And  N_f the number of faces:

    \[N_f = N_b \cdot S \cdot 2\]

Figure 14. How the vertices are placed

To compute the exact position of each vertex, I am using Unity built-in Quaternion.FromToRotation method. Quaternion are used to represent rotations in 3D space. However if I know how to used them, I don’t know about the underlying mathematics, so I won’t detail the vertices creation but instead provide the code behind it:

Vector3[] vertices = new Vector3[(_branches.Count+1) * _radialSubdivisions];
int[] triangles = new int[_branches.Count * _radialSubdivisions * 6];

// construction of the vertices 
for (int i = 0; i < _branches.Count; i++) {
	Branch b = _branches[i];

	// the index position of the vertices
	int vid = _radialSubdivisions*i;
	b._verticesId = vid;

	// quaternion to rotate the vertices along the branch direction
	Quaternion quat = Quaternion.FromToRotation(Vector3.up, b._direction);

	// construction of the vertices 
	for (int s = 0; s < _radialSubdivisions; s++) {
		// radial angle of the vertex
		float alpha = ((float)s/_radialSubdivisions) * Mathf.PI * 2f;

		// radius is hard-coded to 0.1f for now
		Vector3 pos = new Vector3(Mathf.Cos(alpha)* 0.1f, 0, Mathf.Sin(alpha) * 0.1f);
		pos = quat * pos; // rotation
		pos+= b._end; // translation to the end point of the branch
		vertices[vid+s] = pos - transform.position; // from tree object coordinates to [0; 0; 0]

		// if this is the tree root, vertices of the base are added at the end of the array 
		if (b._parent == null) {
			vertices[_branches.Count*_radialSubdivisions+s] = b._start + new Vector3(Mathf.Cos(alpha)* 0.1f, 0, Mathf.Sin(alpha)*0.1f) - transform.position;
		}
	}
}

Faces are then created by connecting the vertices from a branch to its parent. I tried to provide a great schema but it ended up being bad, sorry:

Figure 15. Faces construction, a bad illustration

Let's take a look at our mesh construction.

Figure 16. Mesh construction from the branches list

b. Branches with different size

That's good, but we notice something: all the branches have the same size, while we would like branches to be thin on the extremities and get bigger the closer we get from the root. We want to compute the size of the branches by starting with the extremities. Because extremities were added at the end of the branches list, we can parse the branches from the end to the beginning and compute the size using an invert growth. Let  r_e be the size of an extremity. The size of a branch  s_b can be computed using the following rules:

If the branch have no children, it's an extremity:

    \[    s_b = r_e\]

Else, the size of a branch can be computed as following:

    \[s_b = (\sum_{n=1}^{Nb_{ch}} children[n].{size}^{I_g})^{\frac{1}{I_g}}\]

where  Nb_{ch} is the number of children of the branch
 children[] the array of children of the branch
 I_g inverted growth factor

c. Progressive growth of the branches

Right now, the branches added at the extremities suddenly appear on the screen. We would like our mesh to slowly grow, in a organic way. For that purpose, we will add a bit of code for the computation of the vertices of the extremities. Instead of setting those vertices at the end point of the branch, we will interpolate their position between the starting point and the end point given a factor  t:

    \[  t = \frac{t_{li}}{t_i}\]

where  t_{li} is the time elapsed since the last iteration of the space colonization algorithm
 t_i the time between each iteration

The point  P_e around which we can build the vertices of the extremities can be computed as following:

    \[P_e = S_{branch} + \overrightarrow{D_{branch}} \cdot {Length}\]

where  S_{branch} is the starting point of a branch
 \overrightarrow{D_{branch}} the direction of a branch
 Length the length of a branch

Figure 17. The extremities are now slowly growing between each iteration

Here is another take:

Figure 18. A better render of our growing tree

Source code

The source code for this project can be found on my github: https://github.com/bcrespy/unity-growing-tree

Thank you for reading this article through <3

Follow me on Instagram for more content
References
  1. Modeling Trees with a Space Colonization Algorithm, Adam Runions, Brendan Lane, and Przemyslaw Prusinkiewicz
  2. The Nature of Code, Chapter 8: Fractals, Daniel Shiffman
  3. Spherical coordinate system, Wikipedia

Leave a Reply