How to Render a Fractal, Fast (Part 1)

Recently I've been working on a project I've called Brot, after the Mandelbrot set. It started out as a simple image renderer...

Mandelbrot Set

but soon blossomed into a full-scale video zoom renderer:

The code is in C, and after all my tweaks only comes to 128 lines of code. But before diving into the C, here's a quick introduction to the Mandelbrot set, and how it's calculated. Stay with me—there'll be a little bit of maths, but none of it is too hard. I'm the sort of person who scares easily with complicated formulae, so I'll keep it as simple as possible.

The Mandelbrot Set, as you might have guessed, is a set of numbers. These numbers are complex numbers: they have two components, a real part and an imaginary part. To write these numbers out, we can use this notation: \(a + bi\) This corresponds to a complex number with \(a\) as its real component and \(b\) as its imaginary component. If you're wondering what the \(i\) means, that's simple: \(i\) is defined as the square root of negative one, \(\sqrt{-1}\). This is the simplest imaginary number, and is used to build up complex numbers like we've seen above.

So the Mandelbrot Set consists of complex numbers. For a complex number like \(2 + 3i\) to be in the Mandelbrot set, it has to meet one criterion. Let's see what Wikipedia has to say:

Wikipedia Mandelbrot Explanation

Yikes. That's not too understandable, and it's doing that thing where it uses "i.e." to introduce another, even more confusing way of explaining something. I'll have a go at breaking it down.

For any complex number, we can apply some operations to it and get a new complex number. These operations might be adding, subtracting, multiplying by another number—generally all the things we're familiar with doing to regular numbers, only this time they have two components. If we keep applying these operations, we'll end up with a series of complex numbers. The specific operation we're going to be applying is this:

$$z_{n+1} = {z_{n}}^2 + c$$

with $$ z_{0} = 0. $$

In other words, we start with our complex number to test (\(c\)) and another complex number \(z = 0\). Then, to find the next number in the sequence, we multiply \(z\) by itself, and then add \(c\) to the result. This simple "square-and-add" process is all you need to know to calculate Mandelbrot numbers!

The Wikipedia definition mentions this iteration diverging. In practice, what this means is that no matter how many times we go through the square-and-add process above, the number \(z\) will always stay below some maximum upper limit. If \(z\) stays "bounded" below this upper limit, we can say that our coordinate \(c\) is in the Mandelbrot set.


So the Mandelbrot Set is the set of all complex numbers which stay bounded, however many iterations we apply. But how does this translate into the pretty picture above?

To draw a Mandelbrot Set, we can cheat just a bit—it's hard to draw imaginary numbers on a computer screen. To draw the image, we'll need to divide the screen into pixels, and give each one of those pixels a coordinate. Let's say we want an image at 1080p resolution, or 1920 pixels by 1080 pixels. We give each of these pixels a coordinate from \((0, 0)\) to \((1919, 1079)\), and interpret each coordinate as a complex number: for example, the coordinate \((46, 95)\) might be interpreted by our program as the complex number \(46 + 95i\).

Let's write a C program which, for every \((x, y)\) coordinate in an image, works out whether \((x + yi)\) is part of the Mandelbrot set, and colours that coordinate accordingly. To do this, let's first of all write some boilerplate code to output an image file in .PPM format: it's incredibly easy to work with, and we can convert it to a more modern format with graphicsmagick or similar tools.

And here's the image it produces:

Simple incrementing numbers image

Not too pretty, but hopefully it's clear that get_color keeps incrementing the colour it returns for each pixel, so we get a gradient pattern.

Now, instead of just incrementing that colour, let's use the get_color function to calculate whether the pixel is in the Mandelbrot set or not.

To do this, let's write an iterations function which uses C's inbuilt support for complex numbers to see if our pixel ever goes above some upper limit:

/* Calculate the number of iterations a point takes to leave a bound. */
int iterations(double cr, double ci, int max_it) {  
  double complex z = 0 + 0*I;
  double complex c = cr + ci*I;
  int i = 0;

  for (; i<max_it; i++) {
    z = z*z + c;
    if (cabs(z) > 4)
      break;
  }
  return i;
}

This function takes three arguments: the real component of the complex number cr, the imaginary component ci, and a maximum number of iterations.

After creating two complex numbers c and z, it runs through the algorithm as described earlier, and returns i: the number of iterations before either the point escaped our upper bound, or it hit the maximum allowed number of iterations.

This function is used by our new-and-improved get_color function, which now takes the resolution of the output image as parameters maxX and maxY:

/* Calculate the colour for coordinates (x, y) and write into color[]. */
void get_color(int x, int y, int maxX, int maxY, char color[]) {

  double x1 = linmap(x, 0, maxX, X_LOWER, X_UPPER);
  double y1 = linmap(y, 0, maxY, Y_LOWER, Y_UPPER);
  int its = iterations(x1, y1, MAX_ITERATIONS);

  char r = linmap(its, 0, MAX_ITERATIONS, 0, 255);
  char g = r;
  char b = g;

  color[0] = r;
  color[1] = g;
  color[2] = b;
}

On line 3 of the function, it calls the iterations function with x1 and y1 for the real and imaginary components of the number, and MAX_ITERATIONS as the maximum allowed number of iterations.

N.B.: any time you see capital letters for an identifier, it's probably a constant. In this case, I have a #define MAX_ITERATIONS 255 line at the top of the file. Same goes for X_LOWER, Y_UPPER, etc.

But what are x1 and y1? As it turns out, a range of 0-1919 for the real component and 0-1079 for the imaginary component is a bit rubbish: the mandelbrot set appears absolutely tiny, if it's even visible at all! We'd like to view the classic Mandelbrot image, which looks best when the real component is between \(-2.5\) and \(1.05\), and the imaginary component is between \(-1\) and \(1\).

To do this, I've used the function linmap:

/* Linearly map a value from one range into another. */
double linmap(double val, double lower1, double upper1, double lower2, double upper2) {  
  return ((val - lower1) / (upper1 - lower1)) * (upper2 - lower2) + lower2;
}

linmap is a function I find incredibly useful in almost all maths- or graphics-related projects I work on. It looks complicated, but essentially all it does is convert between scales: if you understand that 40% (40/100) is the same thing as 2/5, then linmap should make intuitive sense:

linmap(40, 0, 100, 0, 5) == 2

The linmap function is also used in calculating our r, g and b values: we want to map the number of iterations (between 0 and MAX_ITERATIONS) to a valid colour (between 0 and 255).

The full, working code is here, but I strongly encourage you to type it all out and play around with the values!

Here's the output of what we just wrote:

Neat! In part 2 we'll cover colouring, and some possible optimizations to make it run even faster.