In parts I and II of the fractal posts, I wrote about the Julia fractal set and its concurrent implementation. In the third post of this series, I am going to talk about its close cousin, the Mandelbrot Set.

Recap - Julia Set

Let us recall the definition of the Julia set. It is a mathematical set of points $z_i$ in the complex plane which satisfies the condition,

$$z_{i+1} = { z_i^2 + c }$$

where $$z_i ≠ ∞$$

We implemented this using the following Python code.

import numpy as np

def julia_set(width, height, zoom=1, niter=256):
    """ A julia set of geometry (width x height) and iterations 'niter' """

    w,h = width, height
    pixels = np.arange(w*h, dtype=np.uint16).reshape(h, w)

    # Pick some defaults for the real and imaginary constants
    # 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

    return pixels

The Mandelbrot Set

The Mandelbrot Set - named after the French mathematician Benoit Mandelbrot who studied its properties extensively, is intimately connected to the Julia set.

It is defined as the set of complex points c, for which the Julia set is connected. In simpler words, this means the set of points where the Julia set contains the origin (0, 0).

The Mandelbrot set then forms the set of co-efficients c of the Julia set, which satisfies this criterion.

Before we implement the Mandelbrot set in Python, let us reimplement the Julia set using complex numbers.

Julia Set - Using complex numbers

Since Python has built-in support for complex numbers, the original Julia Set code can be reimplemented to take advantage of this to be simpler and visually closer to the original mathematical definition of the set. This is the code for the same.

import numpy as np

# Julia Set reimplemented using complex numbers
def julia_set(width, height, cx=-0.7, cy=0.27, zoom=1, niter=256):
    """ A julia set of geometry (width x height) and iterations 'niter' """

    w,h = width, height
    pixels = np.arange(w*h, dtype=np.uint16).reshape(h, w)

    # Pick some defaults for the real and imaginary constants
    # This determines the shape of the Julia set.
    c = complex(cx, cy)

    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)
            z = complex(zx, zy)
            
            for i in range(niter):
                radius_sqr = abs(z)
                # Iterate till the point is outside
                # the circle with radius 2.             
                if radius_sqr > 4: break
                # Calculate new positions
                z = z**2 + c

            color = (i >> 21) + (i >> 10)  + i * 8
            pixels[y,x] = color
  
    return pixels

The changes should be fairly self explanatory to anyone who remembers their complex number algebra and understands the original mathematical definition of the Julia Set.

With this re-implementation, we can now derive the code for the Mandelbrot Set from it.

Mandelbrot Set Generator

The Mandelbrot Set as defined above, is the set of co-efficients c for which the corresponding Julia sets are centered on the origin. The logic of calculating the next point remains the same.

Hence by definition, for every point (x, y) in the real plane, we need to find the set of co-efficients (cx, cy) for the complex co-efficient c in the complex plane which satisfies this criteria.

Let us walk through the Julia Set code step by step and build the corresponding Mandelbrot Set code.

We need to do the imports just as before.

"""
mandelbrot.py - Implement the Mandelbrot set in Python using matplotlib & numpy,
by reusing code from the Julia set.

"""
import matplotlib.pyplot as plt
import numpy as np

The function definition has a slight change as it no longer needs to accept any parameters for the complex constant c, as in this case c always defaults to the origin. However we will use a similar numpy array to hold the pixel data so that part of the code is the same as before.

Note that we also pass two additional arguments namely x_off and y_off for the x and y offsets for the viewport. This is useful when you are zooming the fractal and want to pan to the areas you want to zoom.

def mandelbrot_set(width, height, zoom=1, x_off=0, y_off=0, niter=256):
    """ A mandelbrot set of geometry (width x height) and iterations 'niter' """

    w,h = width, height
    pixels = np.arange(w*h, dtype=np.uint16).reshape(h, w)

Now on to the main logic, the for loop that iterates over the dimensions (w x h).

    # The mandelbrot set represents every complex point "c" for which
    # the Julia set is connected or every julia set that contains
    # the origin (0, 0). Hence we always start with c at the origin.
    
    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
            # We use (x-3*w/4) instead of (x-w/2) to fully visualize
            # the fractal along the x-axis

            zx = 1.5*(x + x_off - 3*w/4)/(0.5*zoom*w) 
            zy = 1.0*(y + y_off - h/2)/(0.5*zoom*h)
            
            z = complex(zx, zy)
            c = complex(0, 0)

The main difference from the Julia Set in this section is that we initialize the co-efficient c to (0, 0) (the origin) every time inside the loop, as per the definition.

Now onto the inner loop that performs the iterations.

            for i in range(niter):
                if abs(c) > 4: break
                # Iterate till the point is outside
                # the circle with radius 2.             
                # Calculate new positions
                c = c**2 + z

In the code for Julia set we iterated over the point z as $ z = z^2 +c$ and stopped when it fell out of the circle with radius 2.

Here since we are interested in the co-efficient space, we need to iterate over that. In other words, the roles of z and c are switched for the Mandelbrot Set, when compared to the Julia Set.

The rest of the code is the same as before, with the difference of the color palette. We want to enhance the blues in this color palette so the code is a bit different.

            color = (i << 21) + (i << 10)  + i * 8
            pixels[y,x] = color
  
    return pixels

Here is the full code, including the function which displays the generated pixels using matplotlib. (Note that we have turned off the display of the co-ordinate axes and also accept a cmap parameter - for the colormap of the plot - in this version of the code.)

import matplotlib.pyplot as plt
import numpy as np

def mandelbrot_set(width, height, zoom=1, x_off=0, y_off=0, niter=256):
    """ A mandelbrot set of geometry (width x height) and iterations 'niter' """

    w,h = width, height
    pixels = np.arange(w*h, dtype=np.uint16).reshape(h, w)

    # The mandelbrot set represents every complex point "c" for which
    # the Julia set is connected or every julia set that contains
    # the origin (0, 0). Hence we always start with c at the origin.

    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
            # We use (x-3*w/4) instead of (x-w/2) to fully visualize
            # the fractal along the x-axis
            
            zx = 1.5*(x + x_off - 3*w/4)/(0.5*zoom*w)
            zy = 1.0*(y + y_off - h/2)/(0.5*zoom*h)
            
            z = complex(zx, zy)
            c = complex(0, 0)
            
            for i in range(niter):
                if abs(c) > 4: break
                # Iterate till the point c is outside
                # the circle with radius 2.             
                # Calculate new positions
                c = c**2 + z

            color = (i << 21) + (i << 10)  + i * 8
            pixels[y,x] = color
  
    return pixels

def display(width=1024, height=768, zoom=1.0, x_off=0, y_off=0, cmap='viridis'):
    """ Display a mandelbrot set of width `width` and height `height` and zoom `zoom`
    and offsets (x_off, y_off) """

    pixels = mandelbrot_set(width, height, zoom=zoom, x_off=x_off, y_off=y_off)
    # Let us turn off the axes
    plt.axis('off')
    # to display the created fractal 
    plt.imshow(pixels, cmap=cmap)
    plt.show()

So much for the code. Let us look at some images.

Fractal Images

Here is the beautiful image of the Mandelbrot fractal that our code generates. The Mandelbrot Set

Let us play around with the code a bit and try and generate different images.

Here is the fractal set zoomed at 2X with an x offset of -300 pixels.

>>> import mandelbrot
>>> mandelbrot.display(zoom=2.0, x_off=-300)

The Mandelbrot Set Fractal Zoomed in at 2X

As you may have expected to see, the fractal is self-similar at higher magnifications.

Here is the fractal with a different color-map than the default one.

>>> import mandelbrot
>>> mandelbrot.display(cmap='magma')

The Mandelbrot Set in magma cmap

You can also play around with the color scheme to generate differently colored images. For example, by changing the color scheme to $ (i<<4) + (i<<12) + i*2 $, you get the following image, with the green color more accentuated around the edges.

The Mandelbrot Set Color Variant

For more information on color scheme experiments with the Mandelbrot Set, visit this discussion.

Concurrent version

It is possible to make this version of the code concurrent and improve its performance by using the multiprocessing module, just like we did for the Julia set in the 2nd part of this series. It is left as an exercise to the reader.

All the code developed in this post is available in our fractals code repository @ github .

Conclusion

I hope you enjoyed this post in the fractal posts series. In future posts, we will explore the mathematical properties of these fractals and also go on to develop code for more interesting fractals.

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