Vanilla Mandelbrot Set renderer in HTML5 canvas and javascript.
This code running in action can be seen at https://csl.name/mandelbrot and you can read about it on https://csl.name/post/mandelbrot-rendering/
Written by Christian Stigen Larsen
https://csl.name
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 function is
C_{n+1} = C_{n}^2 + C_{0}
with the initial condition simply formed by taking the coordinates in the complex plane,
C_{0} = x + iy
Looking at the function, one can easily see that for big initial values, the function should diverge. But for values close to origo (i.e., for |x| and |y| less than 1), we would expect the function to converge to zero, since the product of two numbers less than one will always be less than either of the factors (e.g., 0.5 x 0.4 = 0.2, which is less than both factors).
But if we actually plot it, what we get out isn't any nice plot. Instead, we get an amazingly complex and fractured plot. This is the Mandelbrot set.
You can zoom forever into the plot, and it will present you with an unending complex shape. One can also calculate it's so-called Hausdorff dimension, which yields a noninteger number. Thus, it's a fractal.
Calculating the Mandelbrot set is easy if you do it numerically.
Take any point C_0 = (x, y) and then calculate C_1 = (x+iy)^2 + (x+iy)
and continue doing this. For practical purposes, let's predetermine a
threshold value. If the magnitude of C (defined for complex numbers
as being the distance to origo, or sqrt(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.
If C has not exceeded the threshold value after a predetermined number of
iterations, 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
colors.
The trick is to note that when you calculate Z = Z^2 you'll get values Z, 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:
log(log |Z|) / log 2
This is all explained over at http://linas.org/art-gallery/escape/smooth.html
In my code, I originally used the following smoothing equation:
1 + n - Math.log(Math.log(Math.sqrt(Zr*Zr+Zi*Zi)))/Math.log(2.0);
With some elementary logarithm rules, we can simplify this to
// Some constants
var logBase = 1.0 / Math.log(2.0);
var logHalfBase = Math.log(0.5)*logBase;
// ...
return 5 + n - logHalfBase - Math.log(Math.log(Tr+Ti))*logBase;
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.
We also have to remember that we're using Javascript here, which is a relatively slow language because of its dynamic nature. What's interesting in this regard, though, is to identify performance hot spots in the typical javascript engines. It's interesting to test the code on different browsers.
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 ( sqrt(x^2 + y^2) > threshold ) ...
If we just square the treshold value, we should be able to do away with the square root operation:
threshold_squared = threshold^2
// ...
if ( (x^2 + y^2) > threshold_squared ) ...
You've probably noticed that the plot is reflected vertically over the line
y=0. One can take advantage of this.
The main equation is
C_{n+1} = C_{n}^2 + C_{0}
Setting Cr = Re(C) and Ci = Im(C), we get
C_{n+1} = Cr^2 + 2Cr*Ci*i - Ci*Ci + C_{0}
C_{n+1} = (Cr^2 - Ci^2) + i(2Cr*Ci) + C_{0}
giving us
Re (C_{n+1}) = Cr^2 - Ci^2 + x
Im (C_{n+1}) = 2*Cr*Ci + y
Mag(C_{n+1}) = sqrt(Cr^2 + Ci^2)
If we introduce two variables Tr = Cr^2 and Ti = Ci^2, we get
Re (C_{n+1}) = Tr - Ti + x
Im (C_{n+1}) = 2*Cr*Ci + y
Mag(C_{n+1})^2 = Tr + Ti
Tr = Re(C_{n+1}))^2
Ti = Im(C_{n+1}))^2
So we have now replaced some multiplications with additions, which is normally faster in most CPUs. But, again, this is javascript, and javascript has quite a different performance profile. The code above indeed does not give us any significant speedup --- for a 640x480 image, we only save a hundred milliseconds, or so.
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
offset = 4*x + 4*y*width
so that you can now set the RGBA values as
array[offset+0] = red
array[offset+1] = green
array[offset+2] = blue
array[offset+3] = alpha
There are several ways of optimizing this. For instance, we can simply multiply the whole offset by four, which is the same as shifting all bits left two positions. However, javascript works in mysterious ways, so the customary shift operations may not be as fast as in other languages like C and C++. The reason probably has to do with the fact that javascript only has one data type for numbers, and my guess is that it's some kind of float.
Anyway, we now have
offset = (x + y*width) << 2
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
y*640 = y*512 + y*128 = y*2^9 + y*2^7 = y<<9 + y<<7
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
for ( y=0; y<height; ++y )
for ( x=0; x<width; ++x ) {
// calculate RGBA
var offset = 4*(x + y*with);
image.data[offset + 0] = R;
image.data[offset + 1] = G;
image.data[offset + 2] = B;
image.data[offset + 3] = A;
}
a much faster way would be to do
var offset = 0;
for ( y=0; y<height; ++y )
for ( x=0; x<width; ++x ) {
image.data[offset++] = R;
image.data[offset++] = G;
image.data[offset++] = B;
image.data[offset++] = A;
}
So now we've basically saved the work of doing 2*width*height
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:
- Get an image data array
- For each line: Update the array
- For each line: Copy the array back to the canvas
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:
- Get a line buffer by calling
createImageData(canvas.width, 1) - For each line: Update the line buffer array
- For each line: Call
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
as x, y, steps, threshold), you could in theory calculate all points in
parallel.
Such algorithms are colloquially called embarrassingly parallel.
Now, JavaScript is inherently single-threaded: You can only use so-called green threads, meaning that the javascript engine will yield control between them.
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).
The algorithm is very well suited for vector operations. Most modern computers come with hardware optimizations for such operations (SSE, etc). However, we are again limited to what the javascript engines will do for us.
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.
Copyright 2012, 2018 Christian Stigen Larsen
Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except in compliance with the License. You may obtain a copy of the License at
http://www.apache.org/licenses/LICENSE-2.0
Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.
