Numpy arrays have many powerful features. Some of these are fairly obvious, others are a little more obscure. This article contains demonstrations of how to use a few of these features to assist in plotting fractals, in particular Julia sets and the Mandlebrot set.
Firstly a little background. The Mandlebrot set arises from iterating within the complex plane using the formula:
For the Mandlebrot set Z is initialised at 0, while C corresponds to the co-ordinates under consideration in the complex plane. Points which remain close to 0 indefinitely under iteration are part of the Mandlebrot set, while points which eventually tend towards infinity under iteration are not part of the set. Generally when plotted on computer, the parts of the complex plane which are not included in the Mandlebrot set are colour coded by how many iterations it takes their absolute value to exceed 2 (the point of no return after which they must necessarily keep increasing) giving rise to the attractive images of the set (which is actually the plain coloured object including the origin of the complex plane.
Julia sets cover a wide range of functions . Here we shall just consider those Julia sets which share the same iterative formula as the Mandlebrot set. Julia sets iterate in the same way but are calculated from a different initial condition, that being that Z corresponds to the co-ordinates under consideration in the complex plane while C is a constant (0 for the Julia sets we will consider). Effectively Z and C are swapped for Julia sets compared to the Mandlebrot set. This means that while there is just one Mandlebrot set, there is a different Julia set for every point in the complex plane.
So how do you create a program to plot these fractals. As you may have gathered from the title of this article, I intend to demonstrate some useful features of Numpy which make this nice and easy.
Lets start by considering what packages and variables we will need
'''code to create image files for Julia and Mandlebrot set plots''' import numpy as np # to calculate the images import matplotlib.pyplot as plt # to print the images import os # to save the images import cv2 # to convert the images to movies width, height = 2048, 2048 # image size frames = 1 # 1 produces a single still, higher integers produce a series of images (movies too if > 24) magnification_start = 1024 # zoom factor for first image magnification_end = 1024.0 # zoom factor for final image real_start, imaginary_start = -0.0, 0.0 # centre of coordinates for first image plot real_end, imaginary_end = -0.0, 0.0 # coordinates of final image plot real_seed_start, imaginary_seed_start = -0.5, 0.5 # starting seed for Julia sets real_seed_end, imaginary_seed_end = -0.5, 0.1 # ending seed for Julia sets max_iterations = 200 # the depth of iteration to use for each plot (higher gives more detail but takes longer) # variables controlling set selection and output call_julia, call_mandlebrot = True, False # flags to control which sets are calculated save, display = True, False # flags to control saving and on screeen display options julia_directory = r"C:\Users\Justin\Pictures\JuliaSets\Exp" julia_title = "_julia_set" mandlebrot_directory = r"C:\Users\Justin\Pictures\MandlebrotSets\Exp" mandlebrot_title = "_mandlebrot_set" colour_maps = ['hot', 'cool', 'prism', 'rainbow'] # list of available matplotlib plot colours current_colour_map = 0 # to select from colour_maps
These variable should give us everything we need to create images or even movies of Julia sets or the Mandlebrot set.
Now we need to iterate for each image we want to create, passing the starting conditions into functions to create the images we desire, then we can pass the output of these functions to an image plotting and saving function. Finally we can call another function to stitch movies together if required.
# check if we are processing a still and if so set end values = start values
if frames == 1:
magnification_end, real_end, imaginary_end = magnification_start, real_start, imaginary_start
# cycle through each frame in turn
for i in range (frames):
# set the current scaling for the frame
magnification_current = (magnification_start*(frames-i)/(frames)
+ magnification_end*i/(frames))
# set the current locus for the image
complex_current = complex((real_start*(frames-i)+real_end*i)/frames,
(imaginary_start*(frames-i)+imaginary_end*i)/frames)
# calculate and output the correct types of plot
if call_julia == True:
# create the starting seed for Julia sets
complex_seed = complex((real_seed_start*(frames-i)+real_seed_end*i)/frames,
(imaginary_seed_start*(frames-i)+imaginary_seed_end*i)/frames)
# then change our minds and do the complicated one
a = 2*i*math.pi/frames
complex_seed = 0.7885*math.e**(complex(0, a))
# get our calculated julia plot
output_to_plot = julia_calculate(complex_current,
magnification_current, width, height, max_iterations)
# and plot it
plot_image(output_to_plot, colour_maps[current_colour_map],
display, save, julia_directory, julia_title, i)
if call_mandlebrot == True:
# get our calculated mandlebrot plot
output_to_plot = mandlebrot_calculate(complex_current,
magnification_current, width, height, max_iterations)
# and plot it
plot_image(output_to_plot, colour_maps[current_colour_map],
display, save, mandlebrot_directory, mandlebrot_title, i)
# if we created more than 24 images and saved them we likely want a movie
if frames > 24:
if call_julia == True and save == True:
# make a julia movie
make_movie (julia_directory, julia_title)
if call_mandlebrot == True and save == True:
# make a mandlebrot movie
make_movie (mandlebrot_directory, mandlebrot_title)
OK so now to calculate the functions, this is where Numpy becomes extremely useful. We are going to use a few neat features
- np.linspace to set up some arrays of values representing the real and imaginary values of the points in the complex plane
- np.tile to turn the 1 dimensional real and imaginary arrays into a 2D array of complex values
- create a boolean array to use as a mask (which will itself be updated by masking)
Here is the code for a Mandlebrot function
def mandlebrot_calculate(complex_current, magnification_current, width, height, max_iterations):
'''Calculates a plot of a particular part of the mandlebrot set detemined by the magnification level
and a set of complex coordinates to home in on, the resolution of the plot is determined by width and height
detail is determined by the max_iterations'''
# set up arrays to hold the image data (reshape is needed to avoid (1,) arrays)
x = np.linspace (complex_current.real-(width/magnification_current),
complex_current.real+(width/magnification_current), width).reshape(1, width)
y = np.linspace (complex_current.imag-(height/magnification_current),
complex_current.imag+(height/magnification_current), height).reshape(height, 1)
#initialise Z and C from the Mandlebrot set equation across the matrix
# NOTE That Z and C are almost the reverse of each other for julia and mandlebrot sets
Z = np.zeros((height, width), dtype=complex) # Z starts at 0 for the mandlebrot set (but use complex 0)
C = np.tile(x, (height, 1)) + 1j*np.tile(y, (1, width)) # C is the local complex coordinates for the mandlebrot set
# use a boolean mask to let us know when values have escaped and should not receive furhter updates
mask = np.full ((height, width), True, dtype=bool)
mandlebrot_plot = np.zeros ((height, width))
for i in range (max_iterations):
#calculate Z wherever the mask indicates it has not yet escaped above 2
Z[mask] = Z[mask]*Z[mask]+C[mask]
# update the mask with newly escaped areas above 2
mask[np.abs(Z) > 2] = False
#update the output plot with the latest values
mandlebrot_plot[mask] = i
#increase the contrast between the Mandlebrot set and the last iteration
mandlebrot_plot[mask] = 0
return mandlebrot_plot
While extremely useful, the application of linspace and tile is intuitive and relatively straightforward. Linspace produces an evenly spaced array of values across the interval we specify. We can do this for our real and imaginary ranges and then tile each 1D array on a different dimension to give our values of C for the Mandlebrot calculations. the Z array is even easier to set up as Z starts at zero everywhere for the Mandlebrot set, we just need to ensure that we initialise it as a complex array for future iterations.
The more interesting feature is using boolean masks with Numpy arrays, we shall do this in two places. First we initialise a boolean mask called ‘mask’ setting all values to True. Then when we iterate through the steps of the Mandlebrot calculation. After each iteration we calculate the absolute value of each complex number we have calculated. If it exceeds 2 then we know that point is not part of the Mandlebrot set and can cease calculating, but how can we represent this? Well using a boolean mask achieves this objective rather neatly. The syntax to use a boolean mask is array_to_be_masked[boolean_mask_array]. Rather than conventional ranges being selected as with normal array addressing, the boolean mask refers to all places in the target array where the mask array has a value of True.
First we shall update the mask array by using a mask on it to update escaped points to False:
mask[np.abs(Z) > 2] = False
Once we have updated the mask, we can use it to stencil the iteration number onto our output array using
mandlebrot_plot[mask] = i
This means that each iteration the area exposed by the mask will shrink until eventually it would only contain the points comprising the mandlebrot set (though obviously due to resolution and iteration depth we only get an approximation of this). Then we can pass the completed array of the iteration at which a particular point exceeded an absolute value of two back from the function.
The Julia set function is almost identical, but note that the C and Z set initialisations are reversed. Also note that while the ability to track and zoom is preserved for Julia sets, I have found it preferable to restrict movies to one or other action for Julia sets. By setting the Julia set Seed equal to the focus for a corresponding Mandlebrot set, the interesting visual interrelation between the two sets can be investigated.
def julia_calculate (complex_current, magnification_current, width, height, max_iterations):
'''Calculates a plot of a particular julia set detemined by a particular complex seed
magnification determines zoom, the resolution of the plot is determined by width and heigh
detail is determined by the max_iterations'''
# set up arrays to hold the image data (reshape is needed to avoid (1,) arrays)
x = np.linspace (complex_current.real-(width/magnification_current),
complex_current.real+(width/magnification_current), width).reshape(1, width)
y = np.linspace (complex_current.imag-(height/magnification_current),
complex_current.imag+(height/magnification_current), height).reshape(height, 1)
#initialise Z and C from the Julia set equation across the matrix
# NOTE That Z and C are almost the reverse of each other for julia and mandlebrot sets
Z = np.tile(x, (height, 1)) + 1j*np.tile(y, (1, width)) # Z is the local complex values for Julia sets
C = np.full((height, width), complex_seed) # C is the starting complex seed
# use a boolean mask to let us know when values have escaped and should not receive furhter updates
mask = np.full ((height, width), True, dtype=bool)
julia_plot = np.zeros ((height, width))
for i in range (max_iterations):
#calculate Z wherever the mask indicates it has not yet escaped above 2
Z[mask] = Z[mask]*Z[mask]+C[mask]
# update the mask with newly escaped areas above 2
mask[np.abs(Z) > 2] = False
#update the output plot with the latest values
julia_plot[mask] = i
#increase the contrast between the Julia set and the last iteration
julia_plot[mask] = i
return julia_plot
Now we have the arrays containing the number of iterations before a given point exceeds abs(z)=2. Lets pass these to a plotting function which will use MatPlotLib to generate us some pretty pictures. We don’t need axes or scale bars so this code is reasonably straightforward.
def plot_image(plot_to_show, colour_map, display, save, directory, plot_title, image_number, dpi = 72):
'''Takes an input numpy array and optionally plots (dependent on printout boolean)and saves (dependent on save boolean)
an image via the matplotlib library. colour mappings, save directories and image titles may be set.
The function also takes an image number so that sequential images can easily be saved to make moview'''
fractal = plt.figure()
fractal.set_size_inches(plot_to_show.shape[1]/dpi, plot_to_show.shape[0]/dpi)
axes = plt.Axes(fractal, [0., 0., 1., 1.])
axes.set_axis_off() # we don't actually want to see the axes
fractal.add_axes(axes)
# plot the figure - note we need to flip the picture to have it correct on the Y axis
plt.imshow(np.flipud(plot_to_show), cmap=colour_map)
# save the image if required using serial indexing to allow easy movie making
if save and image_number < 10000:
os.chdir(directory)
if i < 10:
plt.savefig('000'+str(image_number)+plot_title+'.png')
elif i < 100:
plt.savefig('00'+str(image_number)+plot_title+'.png')
elif i < 1000:
plt.savefig('0'+str(image_number)+plot_title+'.png')
else:
plt.savefig(str(image_number)+plot_title+'.png')
if not display:
# hide unwanted plots
plt.close()
This displays or saves our fractal image files as desired. Finally if we produced a lot of image frames then we should make a movie of the resulting images. we shall do this using OpenCV. Note that while OpenCV can be installed via conda install opencv or pip install opencv, the import call in Python is import cv2 even though the latest version is now OpenCV3! Here is a function to create our movies
def make_movie (source_directory, video_title):
'''Given a source directory and a vidoe title, creates an mp4 video comprised of
all the png images in the directory and writes it to the directory'''
os.chdir(source_directory)
video_title = "movie"+video_title+".mp4"
# grab a list of the names of all our images
images = [image for image in os.listdir(source_directory) if image.endswith(".png")]
# take the first image and process into a numpy array
picture_size = cv2.imread(os.path.join(source_directory, images[0]))
# use that to get the dimensions for our video
height, width, channels = picture_size.shape
# Define the codec and create VideoWriter object
codec = cv2.VideoWriter_fourcc(*'mp4v') # Be sure to use lower case
# initialise a new movie
movie = cv2.VideoWriter(video_title, codec, 24, (width,height))
# process each image in turn
for image in images:
# writing it to our movie file
movie.write(cv2.imread(os.path.join(source_directory, image)))
#cv2.destroyAllWindows()
movie.release()
Putting all this together and picking appropriate values for location, magnitude and seed values gives movies like these:
As usual the code can be found on Github. The code is very highly commmented for ease of understanding by those just starting with Python
