24 Apr
2010

Mandelbrot Renderer

By Patrick Ross | No Comments

So a while back, I thought to myself: I’ve always wanted to make a mandelbrot renderer. So I did some research. The basics of how to build a mandelbrot renderer are fairly well explained, but it seems some of the finer points are often missed, so it took me longer than I would’ve liked to build the renderer. Hopefully you can learn from what I discovered.

For a start, we have the basic definition of the mandelbrot set. Graphically, the plane that we render on represents the complex plane, so each point on it is a complex number. Lets call it c. The number c is considered to be part of the mandelbrot set if the complex number zn+1 = zn2 + c, with z0 = 0, remains bounded as n approaches infinity.

This seems fairly straight-forward, however this will only tell us whether a point is in the mandelbrot set or not, so it won’t generate a very nice image for us. It also seems impractical to have to do an infinite number of iterations on the formula to tell whether a point is or isn’t in the mandelbrot set, so we need to firstly introduce some limitations on the algorithm.

The method that we use is called the escape time algorithm. Quite simply, for every point in the region we wish to render, we iterate the formula until we can say for sure that the point is not bounded. Generally this is when the magnitude of z > 2. We use the number, or the “escape time” for the point to give it a colour. We also define a maximum iteration count, say for example 1000 iterations. At this point if the formula still hasn’t “escaped”, we can feel fairly certain that the point will be in the mandelbrot set.

Lets look at some code. If we’re attempting to work out the colour of a point with real component X and imaginary component Y:

// Define the variables
double x = 0;
double y = 0;
int iter = 0;
double xtemp;

We’ll be using x and y as the co-ordinates of z, and xtemp as a temporary variable on each of the iterations.

// Enter the main iteration loop
while ( (x*x + y*y <= 10) && (iter < _max) )
{
 xtemp = x*x - y*y + X;
 y = 2*x*y + Y;
 x = xtemp;
 iter++;
}

Where _max is the maximum number of iterations. These formulae can be fairly easily worked out using complex number math. Note that in my code I define “escaping” as the point where the magnitude of z is larger than sqrt(10).

Now we compare the value of iter to see what colour we should use. If iter == _max, then the point is probably in the mandelbrot set, so it’s coloured black. If iter is any other value, we need to refer to our colouring algorithm.

if (iter == _max)
 return Gdiplus::Color(0,0,0);
else
 return colourOf(iter % _loop); // Discrete colouring

Where _loop is some value that we use to ensure that the colouring loops through a pattern. colourOf is just a function like so:

std::vector<Gdiplus::Color> cols(16);

void initCols(void)
{
 // Sunset gradient
 cols[0] = Gdiplus::Color(25, 7, 26);
 cols[1] = Gdiplus::Color(9, 1, 47);
 cols[2] = Gdiplus::Color(4, 4, 73);
 cols[3] = Gdiplus::Color(0, 7, 100);
 cols[4] = Gdiplus::Color(12, 44, 138);
 cols[5] = Gdiplus::Color(24, 82, 177);
 cols[6] = Gdiplus::Color(57, 125, 209);
 cols[7] = Gdiplus::Color(134, 181, 229);
 cols[8] = Gdiplus::Color(211, 236, 248);
 cols[9] = Gdiplus::Color(241, 233, 191);
 cols[10] = Gdiplus::Color(248, 201, 95);
 cols[11] = Gdiplus::Color(255, 170, 0);
 cols[12] = Gdiplus::Color(204, 128, 0);
 cols[13] = Gdiplus::Color(153, 87, 0);
 cols[14] = Gdiplus::Color(106, 52, 3);
 cols[15] = Gdiplus::Color(66, 30, 15);
}

Gdiplus::Color colourOf(int colour)
{
 return cols[colour];
}

Where the initCols method was called sometime earlier. This is a sunset gradient for those interested. Also note that I’m using Gdiplus colours. The numbers are (in order) R, G, B, so this can fairly easily be applied to any other rendering library.

The only thing left in this rendering process is to convert the x, y coordinates on the screen into X, Y values as their complex number representation. This can be done like so:

double x0 = double(X - _frameWidth / 2) / double(_zoom) + double(_center_x);
double y0 = double(Y - _frameHeight / 2) / double(_zoom) + double(_center_y);

With _frameWidth and _frameHeight as the width and height of the frame to be rendered to, _zoom is the current zoom and _center_x and _center_y are the X and Y values of the center of the rendered frame.

This will generate nice colourful images, however we can make the images BETTER. But how, you ask. Well, the first improvement we can make is to use the smooth colouring algorithm. This algorithm generates smooth gradients on the image, rather than having banded colours. If you haven’t built a copy of the program yet, banded colours look like this:

As you can see, not smooth. The smooth gradient formula that we use replaces the colourOf reference above, and it looks like this.

if (iter == _max)
 return Gdiplus::Color(0,0,0); //255
if (iter == 1)
 return colourOf(1);

int temp = iter;

// Smooth colouring
while (iter < 10)
{
 xtemp = x*x - y*y + X;
 y = 2*x*y + Y;
 x = xtemp;
 iter++;
}
double mag = sqrt(x*x+y*y);
double base = log(log(mag))/log(double(2));
double nsmooth;
if (base < double(iter) + 1)
 nsmooth = double(iter) + 1 - base;
else
 nsmooth = (double)temp;

return colourOfSmooth(nsmooth / 6, _loop); // Smooth linear colouring

In order for the smooth colouring algorithm to give us a reasonable value, generally at least 10 iterations should have been run, possibly more if you want. Dividing by 6 in the last line is just an arbitrary scaling value I used to make the picture look nicer. So what does the colourOfSmooth method look like?

Gdiplus::Color colourOfSmooth(double n, int loop)
{
 int low = (int)n;
 int base = low % loop;
 double offset = n - (double)low;

 Gdiplus::Color lower = colourOf(base);
 Gdiplus::Color upper;
 if (base < 15)
 upper = colourOf(base + 1);
 else
 upper = colourOf(0);

 BYTE R = BYTE(double(upper.GetR() - lower.GetR()) * offset + double(lower.GetR()));
 BYTE G = BYTE(double(upper.GetG() - lower.GetG()) * offset + double(lower.GetG()));
 BYTE B = BYTE(double(upper.GetB() - lower.GetB()) * offset + double(lower.GetB()));

 return Gdiplus::Color(R, G, B);
}

We basically just use smooth the gradient linearly between the two integer values to either side of the smoothed value. Now we have the renderer outputting a nice smooth gradient. There is one further improvement we can make to make our images look better, and that’s adding sampling. By sampling, I mean that for each pixel, we take several “samples” of what the colour would be at different points, and use the average of these as the colour for the pixel. We just use the following algorithm:

double x0 = double(X - _frameWidth / 2) / double(_zoom) + double(_center_x);
double y0 = double(Y - _frameHeight / 2) / double(_zoom) + double(_center_y);
double x1 = double(X + 1 - _frameWidth / 2) / double(_zoom) + double(_center_x);

double sampleSize = (x1 - x0) / sampleWidth;
long R, G, B;
R = G = B = 0;
for (double i = 0; i < sampleWidth; i++)
{
 for (double j = 0; j < sampleWidth; j++)
 {
  Gdiplus::Color temp = _colourOfPoint(x0 + i * sampleSize, y0 + j * sampleSize);
  R += temp.GetR();
  G += temp.GetG();
  B += temp.GetB();
 }
}

return Gdiplus::Color(BYTE(R / sampleWidth / sampleWidth), BYTE(G / sampleWidth / sampleWidth), BYTE(B / sampleWidth / sampleWidth));

Where _colourOfPoint was the method we were using above to get the colour of the pixel.

The entire project can be found in the code section at http://code.space-burger.com/mandelbrot/. If you just want to play around with the program, it is also in the folder.

That’s pretty much it, other than the framework that I used to build the interface.

Here’s some sample high resolution renders created by the program.

–Patrick