The famous Mandelbrot set is a set of points in the complex plane. In essence, what we want to find out is if the iterative function C below will converge to some constant or diverge to infinity.
The definition is simply
zn+1 = zn2 + c
with the initial condition simply formed by taking the coordinates in the complex plane,
z0 = c = x + iy
Pretend for a moment that you've never seen the Mandelbrot plot before. By only looking at the definition above, can you guess how it would look like when rendered? First we need to look at a few expansions of z expressed with x and y.
z0 = x + iy
z1 = (x + iy)2 + (x + iy) = (x2 - y2 + x) + i(2xy + y)
z2 = z12 + (x + iy) = ...
So, as n → ∞, for which values of x and y will zn converge, and for which will it diverge?
Well, for large values of |x| and small for |y| (points close to origo), the sequence might diverge, because x2 would be dominating over y2. Diverged points are painted black, so we can guess that the plot will be black in all directions some distance from origo.
But for points close to origo — say, |x| and |y| less than 1 — we would expect it to converge, because the product of two numbers less than one will always be smaller than each of its parts, giving small values for x2 - y2. So in a circle around origo with a radius of one, we'd expect to see colored pixels.
But what if the signs of x and y differ? Then things quickly get more complicated. In fact, if we plot using a computer, we won't get a nice colored disc centered at origo. What we get is an infinitely complex and fractured plot.
We can even zoom endlessly into the plot and it will still be as non-uniform and complex as before. While a disc would have a dimension of two, a fractal has a so-called so-called Hausdorff dimension which is something in between a line and a plane. We'd get a non-integer dimension; a fractal.
Calculating the Mandelbrot set numerically is easy.
Given the equations above, take any point z0 = (x, y) and then calculate z1 = (x+iy)2 + (x+iy) and continue doing this. For practical purposes, let's decide on a threshold value. If the magnitude of z — the distance to origo, or √(x^2+y^2)— ever becomes larger than this threshold value we will assume that it will diverge into infinity. If so, stop the calculation and plot a black dot at the current location.
For the Mandelbrot set, it can be shown that the threshold value is exactly two, i.e. any sequence with a |zn| > 2 will diverge. Read more about this.
If |z| has not exceeded the threshold value after a predecided number of iterations (which we choose at will), we will assume that the current parameters makes the function converge. In this case, plot a non-black dot at the current location.
I said above that if the function diverges, one should plot a non-black dot. One could simply paint a white dot here. But instead, maybe we want to get an idea of how fast the function is diverging to infinity at this point.
To do this, just take the current value of the number of steps performed and map that against a color spectrum, and paint that color.
So, functions diverging quickly will get about the same color.
If you use the number of iterations to pick a color, you'll get ugly color bands in the plot. There is a really cool trick to get smooth, gradual color changes.
So, you basically calculate
Z = Z^2 until it diverges and make a note of
the iteration count. What we really want, though, is a fractional
iteration count, so we can multiply that with a color value to get smooth
The trick is to note that when you calculate
Z = Z^2 you'll get values
Z^2, Z^4, Z^8 and so on. If you take the logarithm of this, you'll get the
values 1, 2, 4, 8 etc. If you take the logarithm one more time, you'll get
1, 2, 3, 4, 5 etc. So to get a fractional number of iterations, just do:
This is all explained over at linas.org.
In my code, I originally used the following smoothing equation:
With some elementary logarithm rules, we can simplify this to
which is faster. The constant
5 is another little trick, which should
be explained in the code itself.
Finally, when you calculate the color value of a single pixel, it is in reality just the color of a single point in the Mandelbrot set that is situated somewhere inside that pixel.
What I'm saying is that you'll basically get pixel artifacts in the image, especially in dense areas where the color changes (near the black set, for instance).
So what I do is to use random sampling: Just sample a given number of random points inside the pixel and average the sum of the color values. This is equivalent to rendering the plot at a higher resolution and scaling down.
There are many supersampling techniques to use, and the random sampling was chosen because of its simplicity. The problem is that the resulting picture will look a bit blurry (there are ways around this as well).
Calculating the Mandelbrot set is quite slow, but there are a lot of tricks to speed it up.
When speeding up any code, the first step (after making the code correct, of course) is to look at the algorithm and try to use one with a simpler complexity class. Unfortunately, for the Mandelbrot set, this isn't really possible. So the tricks mentioned here are all cases of micro-optimizations. Nevertheless, they will improve the running time quite a lot.
First, let's look at the inner loop. It continually calculates the magnitude of the complex number C, and compares this with a threshold value. Observe that it takes the square root in doing so:
If we just square the treshold value, we should be able to do away with the square root operation:
You've probably noticed that the plot is reflected vertically over the line y = 0. You can take advantage of this mirroring to halve the computation time. I don't, because you'll mostly render plots that are massively zoomed in.
The main equation is
zn+1 = zn2 + c
C = z and
Cr = Re(z) and
Ci = Im(z), we get
If we introduce two variables
Tr = Cr^2 and
Ti = Ci^2, we get
To plot individual pixels in HTML5 canvas, you get an array and you have to calculate the array offset for a given coordinate pair.
I.e., given RGBA pixel format (four positions), an (x, y) coordinate pair and a width and height, you calculate it by
so that you can now set the RGBA values as
Anyway, we now have
Another trick I'd like to mention. Say that the width and height are fixed to, say 640 and 480, respectively. And old trick to multiply y by 640 would be notice that 640 = 512 + 128 = 2^9 + 2^7, giving us
So now we've converted one multiplication into two shifts and an add. In your commodity language and hardware, this might be quite fast in tight innerloops.
Anyway, we still want to be able to use arbitrary heights and widths, so let's skip that one.
By far, the fastest way of accessing the array is by doing it sequentially.
That is, instead of doing
a much faster way would be to do
So now we've basically saved the work of doing
multiplications, or 600 thousand of them, assuming a 640x480 image.
To draw in the canvas, you request an array, update it and copy it back to the canvas.
Of course, you want to reduce the number of such operations. Because we want an animation showing each line as it is drawn, we'll do this:
The trick here, though is to not use
getImageData. You're going to
overwrite all existing image data, so you can use the same buffer for every
line. So instead, we'll use these operations:
putImageData(linebuffer, 0, y_position)to copy only one line
This ensures that we only copy one line per frame update.
Since the algorithm above is referentially transparent, meaning that it
always produces the same output for the same input (where input is defined
x, y, steps, threshold), you could in theory calculate all points in
Such algorithms are colloquially called embarrassingly parallel.
However, there is a new HTML5 APi called web workers that you can use to create real, OS-level threads. That should make it easy to split up plotting into several threads (preferrably one per logical core).
Take a look at the optimizations done to the Mandelbrot set in The Computer Language Benchmarks Game
There are a lot of cool tricks going on there. Most of those use SSE parallelism for hardware speedup or offloads to the GPU.