Fractals can be easily transformed in artistic canvas. In this tutorial I show you how to achieve continuous colouration for the well-known fractals.
A bit of theory A fractal is an object that displays self-similarity on all scales. Mathematically speaking, a fractal is a set generated by a recursive formula in which every new element is calculated starting from the last computed element.
Well known examples of fractals are the Mandelbrot’s and Julia’s sets, the Newton fractal (in the picture above) and the so-called Burning Ship fractal.
Computers are great in repeating the same thing over and over again, so coding a fractal generator program is the best way to explore fractals and their artistic potential. Here are the formulas to generate some notorious fractals:
Mandelbrot’s and Julia’s Sets:
Zn+1 = Zn^2 + cBurning Ship Fractal:
Zn+1 = (|Zr| + i|Zi|)^2 + cNewton Fractal:
Zn+1 = (2 * Zn^3 + 1) / 3Zn^2These formulas teach how to calculate the nth element of the progression, however they don’t explain how to draw a fractal. To do that we need an algorithm.
The Escape Time Algorithm The escape time algorithm is the easiest method to create fractal images. The mechanism is simple: for each pixel of our plot, we iterate the fractal formula and count the loops it takes to the computed value to escape a given threshold. Then we pick a colour according to the number of iterations.
The threshold is composed by an arbitrary number of iterations (the greater, the sharper the image will be) and a fixed constant value equal to 4. You may legitimately ask why exactly 4 and not whatsoever. Obviously, there is a mathematical demonstration for that!
The implementation is straightforward. Here I propose a simple one in C language (in order to keep this snippet accessible to everybody).
// I chose to represent pixels as int[2] and complex numbers as double[2]
// Here are the indexes to access pixel and complex numbers arrays:
# define X 0
# define Y 1
# define R 0 // real part
# define I 1 // imaginary part
c[R] = (pixel[X] - windowWidth) / (zoom * windowWidth);
c[I] = (pixel[Y] - windowHeight) / (zoom * windowHeight);
z[R] = 0.0f;
z[I] = 0.0f;
i = 0;
while (z[R] * z[R] + z[I] * z[I] < 4 && ++i < maxIterations)
{
tmp = z[R] * z[R] - z[I] * z[I] + c[R];
z[I] = z[R] * z[I] + z[R] * z[I] + c[I];
z[R] = tmp;
}
if (i < maxIterations)
{
// get pixel colour
// call to draw function
}
Mandelbrot and Julia sets implementation is basically the same. The only change concerns the initial C constant value. Indeed Julia’s set is a specific case, namely a sub-set, of Mandelbrot’s one. For the Burning Ship fractal we just have to introduce a new concept: the magnitude of a complex number. Its symbol is |Zn| and its value is:
|Zn| = sqrt(z[R] * z[R] + z[I] * z[I])So we have:
// Time Escape Algorithm implementation for Burning Ship Fractal
c[R] = (pixel[X] - windowWidth) / (zoom * windowWidth);
c[I] = (pixel[Y] - windowHeight) / (zoom * windowHeight);
z[R] = c[R];
z[I] = c[I];
i = 0;
while (z[R] * z[R] + z[I] * z[I] < 4 && ++i < maxIterations)
{
tmp = z[I];
z[I] = fabs((double)(z[R] * z[I] + z[R] * z[I] + c[I]));
z[R] = fabs((double)(z[R] * z[R] - tmp * tmp + c[R]));
}
if (i < maxIterations)
{
// get pixel colour
// call to draw function
}Newton fractal is a bit more complicated instead, because we have to deal with powers and derivatives of complex numbers. C language, like C++, Java and others, offers a library to deal with complex numbers. It can be useful for keeping your code cleaner and more expressive. However in this snippet I chose to do without in order to keep the post consistent.
// Time Escape Algorithm implementation for Newton Fractal
z[R] = (pixel[X] - windowWidth) / (zoom * windowWidth);
z[I] = (pixel[Y] - windowHeight) / (zoom * windowHeight);
i = 0;
tmp = 1.0;
while (tmp > 0.000001 && ++i < maxIterations)
{
old[R] = z[R];
old[I] = z[I];
tmp = (z[R] * z[R] + z[I] * z[I]) * (z[R] * z[R] + z[I] * z[I]);
z[R] = (2 * z[R] * tmp + z[R] * z[R] - z[I] * z[I]) / (3.0 * tmp);
z[I] = (2 * z[I] * (tmp - old[R])) / (3.0 * tmp);
tmp = (z[R] - old[R]) * (z[R] - old[R]) + (z[I] - old[I]) * (z[I] - old[I]);
}
// get pixel colour and call to draw functionCountinuos colouring
The real challenge of drawing fractals is to draw beautiful fractals. While sharpening our fractals is quite easy (we just have to increase the iteration threshold), improving their colour smoothness is trickier. There are different techniques to reach a gorgeous colouring effect. The easiest one is to map the escaping values onto a logarithm scale in order to obtain a continuous gradient of iterations.
In simple words: the time escape algorithm returns an integer that tells us how many iterations a pixel needed to escape the loop. An integer is a discrete value. That means that if we want to pass from a value to its next we have to jump. For instance, when we pass from 100 to 101, we skip all the values in between. We ignore 100.1, 100.2 and 100.3 but also 100.33334, 100.33335, 100.33336 and so on.
It is self-evident that if we use a discrete colour scale to compute pixel colours, then our fractal image will have a discrete colouring, that means a lot of choppy and plain colour bands.
To avoid this ugly effect, we can use natural logarithms to transform the magnitude of every escaped pixel into a value between 0 and 1. By adding this value to the iterations count of the considered pixel we obtain at last a continuous scale of iterations values.
continuous_index = iterations + 1 - (log(2) / |Zn|) / log (2)Knowing that a pixel escaped the loop after 105.34 iterations, while another one escaped after 105.38 allows us to pick the colour up in a more precise way.
Endless colour palette repetition
Now that we have a continuous colour index, we just have to do an extra step to find a good colour chart. For our purpose we consider a colour chart to be good when its colours blend smoothly and endlessly (this because fractals can be zoomed indefinitely).
So… Which is the mathematical function that repeats itself continuously and endlessly and that can help us building a perfect fractal colour palette? Sinus and cosinus, of course!
In this wonderful post of Jim Bumgardner you can find all steps to derive the formulas for real time colour-blending. I really suggest you to read in depth the whole post because it is plenty of handy tricks.
channel colour = sin(frequenc * continuous_index + phase) * center + delta)where frequency and phase are arbitrary values, center is the middle of the colour channel range and delta is the maximum variation gap. continuous_index is the value we calculated before. Considering that a whole channel is 255 we obtain:
channel colour = sin(frequency * continuous_index + phase) * 127.5 + 127.5)This is just an example that you can tweak as you prefer. Again, read Jim’s blog carefully! Here are the mixings I used for the pictures above:
// An handy way to operate with RGBA in C language is to represent colours
// as an Union composed by an unsigned integer (the colour itself) and
// an array of 4 unsigned chars (the RGBA channels)
union colour_u
{
unsigned int number;
unsigned char channels[4];
};
union colour_u c;
c.channels[0] = (unsigned char)(sin(0.016 * continuous_index + 4) * 230 + 25);
c.channels[1] = (unsigned char)(sin(0.013 * continuous_index + 2) * 230 + 25);
c.channels[2] = (unsigned char)(sin(0.01 * continuous_index + 1) * 230 + 25);
c.channels[3] = 255; //alpha bit
return c.number;