Visualising the Mandelbrot Set

I suppose it’s about 20 years since I first read about Benoit Mandelbrot and his pioneering work on fractal geometry. Ever since then, I have been intrigued by how fractal images with beautiful and complex textures can emerge from very simple mathematical rules. However, it was not until today that I finally sat down to generate the iconic image of the Mandelbrot Set for myself from scratch In my previous post, I used a short C program to generate a chaotic graph using the logistic map. Mathematically, the logistic map and the Mandelbrot Set are closely related, so it seemed like a good time to give this a try.

The Mandelbrot Set is defined as follows. The following formula can be used to generate a sequence of complex values.

where and are complex values and . A given value of c belongs to the Mandelbrot Set if and only if the sequence remains bounded for all values of n.

I wrote a short program, “mandelbrot.c”, to create an image of the Mandelbrot Set. As it currently stands, the program always generates an image 640×480 pixels in size. The x axis of the generated image is the real part of c and the y axis is the imaginary part of c. Command line arguments are used to specify:

• The real part of the complex value at the centre of the image,
• The complex part of the complex value at the centre of the image,
• The range of real values shown in the image (the range of imaginary values is always scaled proportionately),
• The output image filename (the image is in plain text PGM format).

//
// mandelbrot.c - Visualise the Mandelbrot Set
// written by Ted Burke - last updated 15-12-2012
//
// The command line arguments are the real and imaginary parts of
// the centre value, the real axis range and the output filename.
// 
// To compile with MinGW:
//
//		gcc mandelbrot.c -o mandelbrot.exe
//
// To run (for example):
//
//		mandelbrot.exe 0 0 5 ms.pgm
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define W 640
#define H 480

// Pixel array
int p[H][W];

int main(int argc, char *argv[])
{
	int x, y;
	double real_centre, imag_centre, real_range;
	double cre, cim, cre_min, cre_max, cim_min, cim_max;
	double zre, zim, zre_next, zim_next;
	double limit = 100.0;
	
	// Parse command line arguments
	if (argc != 5)
	{
		printf("Usage: mandelbrot REAL_CENTRE IMAG_CENTRE REAL_RANGE FILENAME\n");
		return 0;
	}
	real_centre = atof(argv[1]);
	imag_centre = atof(argv[2]);
	real_range = atof(argv[3]);
	cre_min = real_centre - 0.5*real_range;
	cre_max = real_centre + 0.5*real_range;
	cim_min = imag_centre - 0.5*H*real_range/W;
	cim_max = imag_centre + 0.5*H*real_range/W;
	
	// Calculate pixel values
	for (y=0;y<H;++y)
	{
		for (x=0;x<W;++x)
		{
			// Calculate c value for this pixel
			cre = cre_min + (cre_max-cre_min)*x/(double)W;
			cim = cim_min + (cim_max-cim_min)*y/(double)H;
			
			// Iterate for current c value
			p[y][x] = 255;
			zre = 0; zim = 0;
			while (fabs(zre) < limit && fabs(zim) < limit && p[y][x] > 0)
			{
				zre_next = zre*zre - zim*zim + cre;
				zim_next = 2*zre*zim + cim;
				zre = zre_next;
				zim = zim_next;
				p[y][x] *= 0.98;
			}
		}
	}

	// Output image to PGM file
	FILE *f = fopen(argv[4], "w");
	fprintf(f, "P2\n# Image\n%d %d\n255\n", W, H);
	for (y=0;y<H;++y)
	{
		for (x=0;x<W;++x)
		{
			fprintf(f, "%03d ", p[y][x]);
		}
		fprintf(f, "\n");
	}
	fclose(f);
	
	return 0;
}

I ran the program with the following command line:

mandelbrot -0.7 0.0 3.5 ms.pgm

which produced the following image:

I then ran the program again with the following command line:

mandelbrot -1.0 0.0 1.0 ms_zoom.pgm

which produced the following image:

To conclude, I wrote a short Python script to generate a sequence of images that I then combined into a single animated gif using ImageMagick’s convert command. Here’s the Python script (filename: “msgif.py”):

#
# msgif.py - Generate an animated Mandelbrot gif
# Written by Ted Burke - last updated 15-12-2012
#

import math
import os

x1, y1, r1 = -0.755625, 0.18125, 4.0
x2, y2, r2 = -0.755625, 0.18125, 0.4
N = 100

rscale = pow(r2/r1, 1/float(N-1))

for n in range(N):
	x = (1 - n/float(N-1))*x1 + (n/float(N-1))*x2
	y = (1 - n/float(N-1))*y1 + (n/float(N-1))*y2
	r = r1 * math.pow(rscale, n)
	
	print "frame %d/%d: " % (n+1, N), x, y, r
	command = "./mandelbrot " + \
		str(x) + " " + str(y) + " " + str(r) + " " + \
		"frame%03d.pgm" % (n)
	os.system(command)

When I ran that script (“python msgif.py”), it produced 100 individual PGM files (“frame000.pgm”, “frame001.pgm”, etc.) which I combined into a single animated gif (“ms.gif”) using ImageMagick’s convert command:

convert -delay 10 frame* ms.gif

The resulting gif animation is shown below:

Advertisements