by Anand | tags: [ fractal julia python tutorial mathematics visualization matplotlib numpy ]
Fractals with Python - Part I
Visualization through code is one of the most creative and beautiful tasks a programmer can do with his craft. In this article and the next couple of articles in this series, we will take a look at creating beautiful fractal images using Python, matplotlib and numpy.
The Julia Set
The Julia fractal set is named after the French mathematician Gaston Julia who investigated its properties during 1915-18.
It generates points on the complex plane for which the series,
$$z_{n+1} = { z_n^2 + c }.$$
does not go to infinity.
Fractal Generation
To generate a fractal image for the Julia Set, we will use the screen to represent the complex plane for the numbers in the set.
We will start with a pixel and use the equation above to iterate the new pixel on the complex plane as long as it is inside the circle with a radius 2 centered around the origin. The color of the pixel will be a function of the number of iterations.
The constant c
will determine the shape of the generated Julia fractal. This constant being in the complex plane, has a real and an imaginary part.
Let us pick the following value for c
.
$$c = (-0.7 + 0.27i)$$
Imagine our starting point p
is (0.5, 1). Here are the calculations.
$$p = (0.5 + i)$$ $$z = p^2 + c$$ $$ = (0.5 + i)^2 + (-0.7 + 0.27i)$$ $$ = 0.25 + 2*0.5i -1 -0.7 + 0.27i$$ $$ = -1.45 + 1.27i$$
The new value of z
is,
$$z = -1.45 + 1.27i$$
So the new distance to origin is,
$$r = \sqrt{1.45^2 + 1.27^2} = 1.93$$
Since this is less than 2, we can iterate to the next z
value. Performing a similar computation as above, you will get the new z as,
$$z = -0.21 - 3.41i$$
The distance to origin is now,
$$r = \sqrt{0.21^2 + 3.41^2} = 3.42$$
The distance is now greater than two and hence the point is outside the Julia set and we stop the iterations. Since we performed 2 iterations for this point, the pixel’s color value becomes 2.
Time to write some code!
The Julia Set Generator
We will now develop our program step by step. First let us get the dependencies out of the way.
We are going to use matplotlib
for plotting the fractal and numpy
for storing the pixel values. Install the following dependencies in your virtual environment.
$ pip install matplotlib numpy
The code we are going to develop is version agnostic so it does not matter whether you are using Python2 or Python3.
Open your favorite editor and create a new file/buffer named as julia.py
. Add the module documentation and the required imports.
"""
julia.py - Implement the Julia set in Python using matplotlib & numpy
"""
import matplotlib.pyplot as plt
import numpy as np
Let us now develop the function which will draw the Julia set. The function will take the width and height of the plot as the main arguments. It will also accept some default keyword arguments namely, the number of iterations
and a zoom
value, defaulting to 256 and 1 respectively.
We will add line numbers to the function code so that it can be referenced easily later.
1 def julia_set(width, height, zoom=1, niter=256):
2 """ A julia set of geometry (width x height) and iterations 'niter' """
3
4 # Function body follows
Next, we create a numpy
array to keep the pixel values. The size of the array should equal the area of the canvas (width*height
) and it should represent a matrix of order height x width
.
5 w,h = width, height
6 # Why (hxw) ? Because numpy creates a matrix as row x column
7 # and height represents the y co-ordinate or rows and
8 # width represents the x co-ordinate or columns.
9 pixels = np.arange(w*h, dtype=np.uint16).reshape(h, w)
Now let us create the constant for the fractal. We will use the same value as before.
10 # The constant "c"
11 # This determines the shape of the Julia set.
12 c_real, c_imag = -0.7, 0.27
Now, we loop over the (x,y) values and transform the cartesian (x,y)
plane to the complex z
plane.
13 for x in range(w):
14 for y in range(h):
15 # calculate the initial real and imaginary part of z,
16 # based on the pixel location and zoom and position values
17 zx = 1.5*(x - w/2)/(0.5*zoom*w)
18 zy = 1.0*(y - h/2)/(0.5*zoom*h)
For each point (zx, zy) in the complex plane, we now iterate using the formula for the set till the point drops out of the circle of radius 2. If we are still within this circle, we iterate to the new point as explained previously.
19 for i in range(niter):
20 radius_sqr = zx*zx + zy*zy
21 # Iterate till the point is outside
22 # the circle with radius 2.
23 if radius_sqr > 4: break
24 # Calculate new positions
25 zy,zx = 2.0*zx*zy + c_imag, zx*zx - zy*zy + c_real
Next, we need to give a color for the pixel. The following formula converts the iteration value to a nice (r,g,b)
value for the pixel. This happens inside the outer for loop.
26 color = (i >> 21) + (i >> 10) + i*8
27 pixels[y,x] = color
Finally, after the loops are done, we display the generated fractal.
28 # display the created fractal
29 plt.imshow(pixels)
30 plt.show()
Calling the function in the Python prompt generates the beautiful fractal as shown below.
>>> import julia
>>> julia.julia_set(1024, 768)
Here is the complete function.
def julia_set(width, height, zoom=1, niter=256):
""" A julia set of geometry (width x height) and iterations 'niter' """
w,h = width, height
# Why (hxw) ? Because numpy creates a matrix as row x column
# and height represents the y co-ordinate or rows and
# width represents the x co-ordinate or columns.
pixels = np.arange(w*h, dtype=np.uint16).reshape(h, w)
pixels = np.arange(w*h, dtype=np.uint16).reshape(h, w)
# The constant "c"
# This determines the shape of the Julia set.
c_real, c_imag = -0.7, 0.27
for x in range(w):
for y in range(h):
# calculate the initial real and imaginary part of z,
# based on the pixel location and zoom and position values
zx = 1.5*(x - w/2)/(0.5*zoom*w)
zy = 1.0*(y - h/2)/(0.5*zoom*h)
for i in range(niter):
radius_sqr = zx*zx + zy*zy
# Iterate till the point is outside
# the circle with radius 2.
if radius_sqr > 4: break
# Calculate new positions
zy,zx = 2.0*zx*zy + c_imag, zx*zx - zy*zy + c_real
color = (i >> 21) + (i >> 10) + i*8
pixels[y,x] = color
# display the created fractal
plt.imshow(pixels)
plt.show()
NOTE: You can find the full code at the repo https://github.com/anvetsu/anviz
Let us now explore the features of this fractal Set.
Exploring the Fractal
We can adjust the various paramters of the Julia set using the function above.
Zooming In/Out
We can zoom in or out of the fractal by passing a different value other than 1
for the zoom
parameter.
>>> import julia
>>> julia.julia_set(1024, 768, zoom=2)
Here is the same zoomed out at 0.5X
>>> import julia
>>> julia.julia_set(1024, 768, zoom=0.5)
Adjusting the Constant
We can change the shape of the fractal by varying the constant. For example, here is the Fractal Set for a value of the constant equal to $$(0.4 - 0.2i)$$
Note how widely this fractal image varies from the original one.
Performance
We can measure the time taken by the code by commenting out the last two lines of the function which performs the plotting and using the timeit
module and building a partial function.
On the laptop that this code was developed - an Octa-core Intel i7 HP Pavilion laptop with 8G RAM running Linux Mint 19 Tara
, the default Julia Set took around 7.7 s, averaged for 10 iterations.
>>> import julia
>>> from functools import partial
>>> julia_default = partial(julia.julia_set, 1024, 768)
>>> timeit.timeit(julia_default, number=10)/10
7.663881219099858
If we increase the size of the fractal, the time taken increases quadratically.
>>> julia_large = partial(julia.julia_set, 1920, 1440)
>>> timeit.timeit(julia_large, number=10)/10
25.005370578500152
The time taken increases with number of iterations, though the growth is lesser than O(n)
.
>>> julia_largeiter = partial(julia.julia_set, 1024, 768, niter=1024)
>>> timeit.timeit(julia_largeiter, number=10)/10
11.184467386000325
Since the fractal is calculated row by row, the performance of our fractal program can be multiplied by many factors by making the calculations concurrent. This can be achieved using different techniques such as multiprocessing
, concurrent futures
and asyncio
in Python.
We can also try and build a fractal explorer where we can adjust the parameters of the fractal and see it getting updated in real-time by using the capabilities of matplotlib
and seaborn
.
Such techniques will be discussed in subsequent iterations of this post. Stay tuned.
References
Note that name and e-mail are required for posting comments