Solarian Programmer

My programming ramblings

Barnsley Fern in Python 3

Posted on November 2, 2017 by Paul

In this article, I will show you how to render the Barnsley Fern in Python 3. The Barnsley Fern is a fractal that can be generated using four simple affine transformations of the form:

\[{w_{i}}(x) = \begin{bmatrix}a_{i} & b_{i} \\c_{i} & d_{i} \end{bmatrix}\cdot \begin{bmatrix}x_{1}\\ x_{2}\end{bmatrix}+\begin{bmatrix}e_{i}\\ f_{i}\end{bmatrix} = {A_{i}} \cdot x + {t_{i}}\]

where the coefficients of the transform are:

     w           a           b           c           d           e           f           p     
1 0 0 0 0.16 0 0 0.01
2 0.85 0.04 -0.04 0.85 0 1.6 0.85
3 0.2 -0.26 0.23 0.22 0 1.6 0.07
4 -0.15 0.28 0.26 0.24 0 0.44 0.07

In the above table p represents the probability factor for a transform. For example, the second transform will be used 85% of times, third transform 7% and so on.

From a practical point of view, the Barnsley Fern is generated starting with an initial point \({x_{1}}\), \((0, 0)\), and iteratively calculating the next point \({x_{i}}\) using one of the above \({w_{i}}\) transforms.

At each step, we generate a random number r from the interval [0, 1) and, interpreting this number as a probability, we pick the corresponding transform. For example, if r is in the interval \([0, 0.01]\) we pick \({w_{1}}\), if r is in the interval \((0.01 + 0.85, 0.01 + 0.85 + 0.07]\) we pick \({w_{3}}\) and so on. Note, we define the intervals using the cumulative sum of the probability factors.

This is the result, after 1000000 iterations:

Barnsley Fern

Now, we can implement the Python code that will generate the above fractal. Because we can generate an infinite number of fractals, by varying the coefficients of the affine transforms and their probability factor, I propose to implement a general Barnsley class that will let us use more than a set of coefficients. You can find the entire code on my GitHub https://github.com/sol-prog/Barnsley-Fern.

 1 class Barnsley:
 2     """Generate Barnsley like fractals using affine transforms"""
 3     def __init__(self, nr_points, coefficients="fern"):
 4         self.nr_points = nr_points
 5 
 6         #Initial (starting point)
 7         self.x = 0.0
 8         self.y = 0.0
 9 
10         # Store the fractal points
11         self.point_x = []
12         self.point_y = []
13 
14         # Select the set of coefficients to use
15         if coefficients == "fern":
16             self.probability_factors = [0.01, 0.85, 0.07, 0.07]
17             self.a = [0, 0.85, 0.20, -0.15]
18             self.b = [0, 0.04, -0.26, 0.28]
19             self.c = [0, -0.04, 0.23, 0.26]
20             self.d = [0.16, 0.85, 0.22, 0.24]
21             self.e = [0, 0, 0, 0]
22             self.f = [0, 1.6, 1.6, 0.44]
23 
24         # other sets of coefficients are present in the code from GitHub:
25         # https://github.com/sol-prog/Barnsley-Fern
26 
27 
28         self.nr_transforms = len(self.probability_factors)
29 
30         # Cumulative sum of the probabilty factors,
31         # this defines the intervals corresponding to each transform
32         self.cumulative_probabilities = [0] * (self.nr_transforms + 1)
33         for i in range(1, len(self.cumulative_probabilities)):
34             self.cumulative_probabilities[i] = self.cumulative_probabilities[i - 1] + \
35                                                self.probability_factors[i - 1]

Next, we need a function to select which transform to use:

1     def select_transform(self):
2         """Randomly select an affine transform"""
3         rnd = random.random()
4         for i in range(self.nr_transforms):
5             if self.cumulative_probabilities[i] <= rnd <= self.cumulative_probabilities[i + 1]:
6                 self.current_transform = i
7                 break

We can calculate and store the next point of the fractal using:

1     def next_point(self):
2         """Get the next point of the fractal"""
3         self.select_transform()
4         x_new = self.a[self.current_transform] * self.x + self.b[self.current_transform] * self.y + self.e[self.current_transform]
5         y_new = self.c[self.current_transform] * self.x + self.d[self.current_transform] * self.y + self.f[self.current_transform]
6         self.x = x_new
7         self.y = y_new
8         self.point_x.append(x_new)
9         self.point_y.append(y_new)

Use the above function to iteratively generate the number of points you need for the fractal:

 1     def generate_points(self):
 2         """Generate all the fractal points"""
 3         for _ in range(self.nr_points):
 4             self.next_point()
 5 
 6         # Bounding box for the fractal
 7         self.x_min = min(self.point_x)
 8         self.x_max = max(self.point_x)
 9         self.y_min = min(self.point_y)
10         self.y_max = max(self.point_y)

Next step is to visualize the generated fractal. A simple solution is to use an image file for this.

In my last article I’ve shown you how to write your own PPM image generator from scratch in Python and you can follow this approach if you want. In the present article I will show you how to use the Pillow image library. On most systems you can install Pillow with:

1 pip install pillow

Next, we initialize the fractal points, transform these points to the image space and save the image:

 1 def main():
 2     """Generate and the save as a PNG image a Barnsley Fern"""
 3 
 4     # Initialize the fractal data
 5     nr_points = 100000
 6     fern = Barnsley(nr_points)
 7     fern.generate_points()
 8 
 9     # Define the image size and scale factor for the fractal data
10     width, height = 500, 500
11     scale = min([height/(fern.y_max - fern.y_min), width/(fern.x_max - fern.x_min)]) * 0.9
12 
13     # Initialize an array that will store the image pixel data
14     image_data = array.array('B', [255, 255, 255] * width * height)
15 
16     # For every point of the fractal data, transform the point in the image space
17     # and fill the pixel color
18     for i in range(nr_points):
19         x = int((fern.point_x[i] - fern.x_min) * scale) + int((width - (fern.x_max - fern.x_min) * scale)/2)
20         y = -int((fern.point_y[i] - fern.y_min) * scale) - int((height - (fern.y_max - fern.y_min) * scale)/2)
21 
22         index = 3 * (y * width + x)
23         image_data[index] = 0
24         image_data[index + 1] = 255
25         image_data[index + 2] = 0
26 
27     # Show and save the image
28     img = Image.frombytes("RGB", (width, height), image_data.tobytes())
29     img.show("Barnsley's Fern")
30     img.save("barnsley_fern.png")

As previously mentioned, you can generate an infinite number or fractals by altering the coefficients of the transform. Here is what you get with a slightly changed set of coefficients:

     w           a           b           c           d           e           f           p     
1 0 0 0 0.16 0 0 0.04
2 0.7 0.035 -0.04 0.8 0 1.6 0.8
3 0.2 -0.29 0.23 0.22 0 1.6 0.08
4 -0.2 0.28 0.26 0.25 0 0.44 0.08

Custom Fern

If you want to read more about the mathematics of fractal generation, I recommend reading Fractals Everywhere by Michael F. Barnsley:

If you want to learn more about Python, I recommend reading Python Crash Course by Eric Matthes:


Show Comments