Intro to fractal art with Python

A low-level introduction to drawing fractals on the complex plane

11 min readFeb 26, 2022

If you own a computer, drawing fractals is simple (otherwise not so much). In this post, we’ll see how two simply-implemented algorithms allow you to draw an unlimited number of different fractals on the complex plane. Throughout all this, only a bare minimum knowledge of mathematics and programming is needed (if you know how distances work, what complex numbers are, and that multiplying and adding them results in a polynomial, you’re good to go).

The code featured in this post can be found here. It also includes a script called fractal.py that implements the ideas of this post into a CLI app.

What you’ll need

For this, you will only need python installed on your computer, along with the libraries matplotlib, NumPy and OpenCV. And maybe familiarize yourself with the RGB(A) color format.

What are fractals?

First, we need to get an idea of what we’ll be drawing. What even is a fractal? We often say that something is a fractal if its smaller parts resemble the whole part.

But here we’ll focus on two types of specific fractals: Julia sets and the Mandelbrot set. They are all defined as subsets of the complex plane.

Julia sets

The Julia set is defined for any given complex polynomial p, as the set of complex numbers z, for which it holds that by applying p iteratively, the process stays bounded. This means that the sequence

z, p(z), p(p(z)), p(p(p(z))), …

can be fitted inside a large enough circle on the complex plane.

To decide whether a point is inside the Julia set, we’ll need to determine if the above sequence stays bounded. For this, we have to use a rather technical theoretic result.

It tells us there exists such a number r that for all complex numbers z outside the circle with radius r (so |z| > r), the above sequence will be unbounded. This also implies that if our sequence at any point escapes the circle, iteration will be unbounded from then on.

The Mandelbrot set

First, we have to define the following polynomial:

q(c, z) = z² + c

It can be viewed as a family of polynomials, each corresponding to its own complex number c. The Mandelbrot set is the set of all complex numbers c, for which the polynomial q(c, _) iterated from 0 stays bounded. This means

0, q(c,0), q(c, q(c,0)), q(c, q(c, q(c,0)), …

can be fitted inside some complex-plane circle.

The (non-trivial) motivation behind the Mandelbrot set is that it is a family of all the points c, for which the polynomial q(c, _) has a connected Julia set.

Similarly, as with Julia sets, we want to test whether a given c produces an unbounded sequence. It can be shown that in this case, we can take r=2.

Escape-time

So, to draw fractals we need to decide whether iterated polynomial sequences remain bounded. Escape-time is a simple algorithm that tests this by generating the iterative sequence and counting the steps needed for a term to escape outside a given radius R (should be set to what was previously denoted r).

We represent polynomials with lists and evaluate them using the horner algorithm

If the algorithm returns k < K, we know the sequence escaped our radius. This is proof that the sequence is unbounded. But what if we get k = K? Then the sequence remained inside the circle for all of the K iterations. If K was large enough, we have a right to assume the sequence is bounded!

Drawing the Mandelbrot set by separating those two options (and setting R=2) would look like:

Ew, that was pretty gross. Where are all the colors? We have to see that the algorithm does more than what we’ve used it for. It tells us the TIME it took for our sequence to ESCAPE the circle of radius R. We use that to draw the set a bit more colorfully. But first, let’s meet some tools in python that make these things easier.

Matplotlib’s colormaps

Let’s focus on what we have so far. We have an algorithm that assigns a number from 1 to K to a given point z in a complex plane. So if we map those numbers to colors, we get a mapping from the complex plane to the set of colors. In other words, we color the complex plane and get a picture.

We can define any color map we want, but the easiest way is to use matplotlib, which comes with predefined color maps (or lets you define your own). Though note that they map from the interval [0,1] to [0,1]⁴ (the codomain are tuples (red, green, blue, alpha)). The way to fix this is to norm our colors, so instead of i, we pass to the function i / K.

OpenCV image format

Images will be represented with NumPy arrays, where in the (i,j) position, we’ll encode the color of the (i,j)-th pixel in BGR (reversed RGB) format. For example, some 5x5 pixel image would look like this:

155 15 1  155 15 1  169 22 1  169 22 1  169 22 1  
155 15 1  183 29 2  183 29 2  206 45 4  234 78 13  
155 15 1  195 37 3  235 211 57  0 0 0  0 0 0  
155 15 1  195 37 3  235 211 57  0 0 0  0 0 0  
155 15 1  183 29 2  183 29 2  206 45 4  234 78 13

To make our life easier, we can define a function that applies a given mapping to each coordinate, assigning it a BGR value, and then draws the image using the OpenCv package.

Drawing the Mandelbrot set

All that’s left is to map our (i,j) pixel coordinates to values on the complex plane. We simply do this by slicing a complex-plane square (with a given edge-length 2r and centre ctr) into an nxn grid, via the function mapToComplexPlaneCenter.

We used 2000² pixels and K=80 allowed iterations to produce the below image.

Still ew for many reasons. What are all those lines doing there? It’s supposed to be smooth. The problem is that our algorithm returns natural numbers, which means its codomain is a bit rough. We’ll solve this later on.

Drawing Julia sets

Let’s pick some polynomial, say

p(z) = z² + (0.285 + 0.01i)

Remember that we’ve said there exists some r, telling us when we’ve iterated into oblivion. Its definition is in the code (as radiusJulia).

The lines we are trying to get rid of can often be pretty awesome though (if you define your color map appropriately you can make anything look awesome).

Am beating myself up over losing the polynomial that produced this.

Generally, don’t expect to find an amazing Julia set on any random polynomial. They can be hard to find. Some polynomials that produce interesting Julia sets can be found online or through an explorer.

Distance Estimation Method (DEM)

Let’s try to solve the aesthetic problems we’ve encountered before. This will be done by implementing DEM. Whereas escape-time depends on somewhat trivial results, we won’t pretend to understand DEM, as it depends on some highly non-trivial theorems (e.g. the Koebe 1/4-theorem, Douady-Hubbard theorem). DEM estimates the distance between a given point and a fractal set. This improves escape-time, as distance as a function is continuous and thus produces nice gradients without any visible lines.

DEM for the Mandelbrot set

The algorithm goes a bit like this:

A slight issue in applying DEM is that we don’t know what the largest number returned can be and thus can’t norm our result to fall into [0,1], where our color mapping operates. This is solved by saving all of our values and then normalizing them as a whole (dividing them by their maximum element).

We optimized DEM a bit and had to tweak some things in order to get better pictures. E.g. we took a minus logarithm of the distance returned by DEM. This is done to achieve a more noticeable color contrast near the set boundary.

We can also zoom closer into some specific parts (look at the optimization chapter for recommendations on the places to zoom into).

I don’t remember where this was but I used a reversed `gist_stern` colormap.

When you zoom in on a set, remember to increase the iteration count K, as being closer to the set’s boundary requires better approximations.

DEM for the Julia set

Here, the algorithm is fairly similar.

The only difference is we also require an argument dp, which is the derivative of polynomial p. It is easy to compute, as seen in the below implementation.

Our code also applies a function over the whole normalized array of values returned by DEM (in the above code, we use pow on the array normalized to get the array adjusted).

This is useful when many colors are hidden close to the border of the set, so the normalized numbers are either very near 1 or very near 0. We can use the fact that applying a (positive) power of less than 1 brings numbers in (0,1) closer to 1 and applying a power of more than 1 brings them closer to 0.

In the above case, we’ve slightly brought purple out of hiding near the boundary by applying a power of 0.8. Compare with the image generated without using this technique:

It’s especially useful, when we have a really thin set and the only visible points are those indicating the set to be close. In that case, we can push some percent of the highest numbers even higher (closer to 1) and the rest closer to 0. This is easily done.

Using this on the polynomial

z² + (-0.7510894579318156 + 0.11771693494277351i),

with color map cubehelix and center -0.5 with radius around 0.5, we get something like the image below

This example (not zoomed) is also shown in the repo’s readme.

Optimizations

As so far we’ve prioritized readability, the code was very inefficient. There are many things to consider when optimizing. If you can optimize a function that runs for every pixel, it shows. For instance, we make many function calls to modularize our code, but they are unnecessary and cost us time. We also run one process, whereas the code is easy to parallelize.

The big issue of course is Python, which is slow and used only for its readability and easy image generation. But Python is implemented in a much faster language - C and we can write the heavy iteration part of the code in C and import it in Python (that’s for another post).

A more theoretic optimization is realizing that if a point lies within a specific cardioid or circle, then it’s in the Mandelbrot set. Specifically, c has to satisfy

We can avoid the square-rooting and optimize the checking of those two conditions.

An interesting fact is that the points that satisfy equality in the above equation lie at the boundary of the Mandelbrot set. So those points may be a good choice to zoom into.

Cardioid and circle inside the Mandelbrot set

Even more, if we remember what the Mandelbrot set represents, we see that at the border of the set are such values c, for which the polynomial p(z) = z² + c has a close-to unconnected Julia set. So, polynomials of that form may produce interesting Julia sets.

Exercise: find the c for which **p(z) = z² + c** produces this set

Further applications

Some ideas for further implementations:

Coloring the fractal’s interior. Throughout the post, we kept the color of the fractal interior constant. We can change that, e.g. by color-mapping the distance made by a complex number during the iteration.

Different representations of a complex plane. Our mapToComplexPlaneCenter was an affine mapping to a complex plane, which meant it uniformly sliced the square of the complex plane, where we wanted to draw. We could make it less natural and e.g. use squaring.
Assigning fractal images to objects. We know how to map color maps and polynomials (or points with radiuses) to fractal images. Now we can map any other object to a color map and a polynomial so that by composition we get a mapping from that object to a fractal image. That object can be e.g. a name, a file, …
Make a fractal video. We can parametrize a family of polynomials, indexed by t and view t as time, going through some interval [a,b]. Here is an example.
Plot 3d fractals. This is way harder for many reasons. Also, complex numbers have no natural 3d generalization. But it is possible to mimic complex operations in 3d. That’s the idea with the Mandelbulb fractal. Another idea is to use quaternions, which are a 4d generalization of complex numbers. Overall, the best site I can refer you to is this one.

Happy fractaling!

Some actual math

For mathematically inclined, here we state and prove the theorems making escape-time valid, by defining the escape-radiuses.

We iteratively apply the proposition and see that our Julia-sequence becomes unbounded. Note that the above max{...} is what we’ve denoted by r at the beginning and named radiusJulia in the code.

And that’s it. Just so we don’t end all technical, we also state perhaps the most beautiful theorem of this theory. It’s called The Fundamental Dichotomy and says that every Julia set is either connected or has uncountably-many different connected components (if you don’t know, uncountable means a lot).

Conclusion

We’ve seen how to draw acceptably pretty fractals using not-really-that-complicated tools. It’s awesome that just simple iterations produce such complexity and take symmetry to another level. Let’s finish symmetrically.