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)

The Julia Set

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)

The Julia Set Fractal Zoomed in at 2X

Here is the same zoomed out at 0.5X

>>> import julia
>>> julia.julia_set(1024, 768, zoom=0.5)

The Julia Set Fractal Zoomed out at 0.5X

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)$$

The Julia Set Fractal with a different c

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

  1. http://paulbourke.net/fractals/juliaset/
  2. https://lodev.org/cgtutor/juliamandelbrot.html

Note that name and e-mail are required for posting comments