profile picture
07 May 2015
Fractals continuous colouring

An efficient and smooth coloration tecnique is the best way to trasform computer generated fractals in artistic canvas. In this short tutorial I'm going to introduce an easy way to achieve a continuous coloration of the most well-known fractals.

Newton Fractal
Newton Fractal

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 + c
Burning Ship Fractal:
Zn+1 = (|Zr| + i|Zi|)^2 + c
Newton Fractal:
Zn+1 = (2 * Zn^3 + 1) / 3Zn^2

These 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 color 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 numebers 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 color
  // call to draw function
Mandelbrot's set
Mandelbrot's Set

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 color
  // 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 color and call to draw function

Countinuos coloring

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 color 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 color scale to compute pixel colors, then our fractal image will have a discrete coloring, that means a lot of choppy and plain color 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 color up in a more precise way.

Julias's set
Julia's Set

Endless color palette repetition

Now that we have a continuous color index, we just have to do an extra step to find a good color chart. For our purpose we consider a color chart to be good when its colors 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 color 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 color-blending. I really suggest you to read in depth the whole post because it is plenty of handy tricks.

channel color = sin(frequence * continuous_index + phase) * center + delta)
where frequence and phase are arbitrary values, center is the middle of the color 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 color = sin(frequence * continuous_index + phase) * 127.5 + 127.5)
This is just an exemple 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 colors
// as an Union composed by an unsigned integer (the color itself) and 
// an array of 4 unsigned chars (the RGBA channels)

union       color_u
  unsigned int  number;
  unsigned char channels[4];

union color_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;