How to plot and visualize Mandelbrot set and Mandelbulbs in Datalore

Datalore's built-in tools are great for interactive plots and charts - but our visualization libraries are capable of so much more. With just a handful of Python lines, you can create captivating images of fractals.

Visualizations of the set of complex numbers known as the Mandelbrot set result in intricate fractal-like images which drew attention to the Mandelbrot set outside of mathematics. Inspired by Syntopia's series of posts, today we're plotting the Mandelbrot set in 2D and 3D using the datalore.plot library (which implements the well-known R library ggplot).

You can play with zooming and constructing fractal sets in this Datalore workbook. All Mandelbrot visualizations in Datalore start with a single function - plot_mandelbrot_set.

Mandelbrot set

Let's start with the definition. Consider a complex number $c\in\mathbb{C}$ and the following recursive formula: $\displaystyle z_{n+1}=z_n^2 + c$, where $z_0=0$. If $|z_n|<\infty$ as $n\rightarrow\infty$, then $c$ belongs to a Mandelbrot set $\mathcal{M}$.
Note that you can consider $z_n$ as a function $f(n, c)$ such that $f(n+1, c) = f^2(n, c) + c$. In practice we'll assign point $c$ to $\mathcal{M}$ if $|z_n|>H$, where $H$ is some arbitrary value (greater than 2) called horizon value.
One can get a monochrome picture of the Mandelbrot set by just checking whether or not the points belong to $\mathcal{M}$, but you might remember colorful zoomings of this famous fractal set.
One of the ways to assigns color to an exterior point while plotting a Mandelbrot set is to use escape time values - the number of steps needed for a length $r_n=|z_n|$ to achieve horizon level $H$.

We can perform the assignment process with the help of the following function:

def mandelbrot(c, threshold=2, num_iter=32):
    z = c
    for i in range(num_iter):
        if np.abs(z) > threshold:
            return i
        z = z*z + c
    return 0

Here we assume that $c$ is a complex number:

c = complex(x, y)

Laplace operator

Let's plot a monochrome Mandelbrot set and the result of applying a Laplace operator to an image $u(x,y)$ of the Mandelbrot set.

Laplace operator:
$$\displaystyle\Delta u=\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}$$

Finite-difference version of Laplace operator:
$$\displaystyle\Delta u\simeq \frac{u(x+h,y)+u(x-h,y)+u(x,y-h)+u(x,y+h)-4u(x,y)}{h^2}$$

first

Anti-aliasing

We can see that the color changes stepwise on our picture. To smooth the transition between color segments, we will apply anti-aliasing - a renormalization technique for image rendering.

Note how that integer-valued escape time in the above function depends on the horizon value (threshold argument of our function). Our goal is to get real-valued escape time which will not depend on actual horizon value.
Define the function $\displaystyle f(x)=\frac{\log\log(x)}{\log 2}$ and let $t$ be the escape time and $H$ the horizon value. Then the renormalization we are looking for takes the form:
$$t - f(|z|) + f(H)$$
With this renormalization, we obtain the following Mandelbrot image.

second

Histrogram equilization and adding color

You might have noticed that we have a few step-wise color changes close to the Mandelbrot set itself. One way to mitigate this problem is to perform color equalization with the help of the power function:
$$\left(t - f(|z|) + f(H)\right)^w$$
for some $w < 1$. The following function was used to produce the image of the Mandelbrot set:

def mandelbrot_set(x_min, x_max, y_min, y_max, width=270, height=270, threshold=2**40, num_iter=30, power=0.2):
    xs = np.linspace(x_min, x_max, width)
    ys = np.linspace(y_min, y_max, height)
    img = np.array([mandelbrot(complex(x, y), threshold, num_iter)**power for y in ys for x in xs])
    return img

To add colors, assign color intensities to the RGB channels of your image:

def plot_colored_image(image, rgb, shape, x_size, y_size, title):
    width, height = shape
    colored_image = np.zeros((width, height, 3))
    red, green, blue = rgb
    colored_image[:,:, 0] = image * red
    colored_image[:,:, 1] = image * green
    colored_image[:,:, 2] = image * blue
    plot = ggplot() + geom_image(colored_image) + ggsize(x_size, y_size) + ggtitle(title)
    return plot

third

Zoom in Mandelbrot set

A simple helper function will recalculate the image frame for the selected center and zoom of the picture:

def get_boundaries(center, span, zoom):
    return center - span/2.**zoom, center + span/2.**zoom 

fourth
fifth
sixth

Mandelbulb

There is no strict analogy for complex numbers in 3D space, so there are a number of ways to determine 3D fractals. Some of them use quaternions.

Consider 3D polar coordinates:
$
$\begin{cases}
r = \sqrt{x^2+y^2+z^2}\\\theta = \arctan\left(\frac{\sqrt{x^2+y^2}}{z}\right)\\\phi=\arctan\left(\frac{y}{x}\right)
\end{cases}$
$
that transfer to the Descartes coordinates in the following manner:
$
$\begin{cases}
x = r\sin\theta\cos\phi\\y = r\sin\theta\sin\phi\\z = r\cos\theta
\end{cases}$
$

Given $v = (x,y,z)$, we can define the $n$-th power of this triple:
$$(x,y,z)^n=r^n\left(\sin n\theta\cos n\phi, \sin n\theta \sin n \phi, \cos n\theta\right)$$
Now we can define the recurrent formula for the 3D Mandelbrot set (called a Mandelbulb) as $\displaystyle v_{n+1}=v_n^p + c$, where $c=(x,y,z)$. A typical choice for the power $p$ is 8.

Projection plane

Our goal is to plot a Mandelbulb. Suppose the observer position is defined by a point $Q=(a,b,c)\in\mathbb{R}^3$.
Consider a plane $\Omega$, passing through the origin with a normal vector passing through $Q$: $\displaystyle a\cdot x + b\cdot y + c \cdot z=0$
Suppose $c\neq 0$, then fixing $x$ and $y$ one can derive $z$: $\displaystyle z =-\frac{ax+by}{c}$
The following function will return points belonging to the plane $\Omega$:

def get_plane_points(Q, center, span, zoom, width, height, eps=1e-4):
    x_min, x_max = get_boundaries(center[0], span[0], zoom) 
    y_min, y_max = get_boundaries(center[1], span[1], zoom)
    a, b , c = Q
    x = np.linspace(x_min, x_max, width)
    y = np.linspace(y_min, y_max, height)
    x, y = np.meshgrid(x, y)
    x, y = x.reshape(-1), y.reshape(-1)
    if np.abs(c) > eps:
        z = -(a*x + b*y)/c
        P = np.vstack((x, y, z)).T
    elif np.abs(a) > eps:
        z = -(c*x + b*y)/a
        P = np.vstack((z, y, x)).T
    elif np.abs(b) > eps:
        z = -(a*x + c*y)/b
        P = np.vstack((x, z, y)).T
    return P

Ray tracing algorithm

Starting from the point $Q$ we are going to move along rays in the direction of points $P\in\Omega$. In each new position, we need to know the distance to the closest point belonging to the Mandelbulb $\mathcal{M}$.

Suppose we know how to estimate this distance $r$. This means that we can move one more step in this direction with a step of size $r$. We stop moving in a particular direction if the current distance estimate is below the minimal distance defined by a user.

Here's the function for such an algorithm:

def trace(start, directions, max_steps, min_distance, iterations, degree, bailout, power):
    total_distance = np.zeros(directions.shape[0])
    keep_iterations = np.ones_like(total_distance)
    steps = np.zeros_like(total_distance)
    for _ in range(max_steps):
        positions = start[np.newaxis, :] + total_distance[:, np.newaxis] * directions
        distance = DistanceEstimator(positions, iterations, degree, bailout)
        keep_iterations[distance < min_distance] = 0
        total_distance += distance * keep_iterations
        steps += keep_iterations
    return 1 - (steps/max_steps)**power

Direction vectors

Consider a point $P\in\Omega$ and a start point $Q$. Define direction vector as a unit vector pointing from $Q$ to $P$.
Note that we can vectorize this function:

def get_directions(P, Q):
    v = np.array(P - Q)
    v = v/np.linalg.norm(v, axis=1)[:, np.newaxis]
    return v

Distance estimator

It was shown that the distance to the closest point can be estimated via $$r = \lim_{n\rightarrow\infty}\frac{|v_n|\log |v_n|}{|v_n^{\prime}|}$$
This means that we have to keep track of distances and their derivatives in our algorithm. Denote $r_n=|v_n|$ and $dr_n=|v_n^{\prime}|$.
The most wonderful thing is that we can just use the scalar version of the formula for the derivative update:
$$dr_n=p\cdot r_{n-1}^{p-1}dr_{n-1} + 1$$
The following function summarizes this algorithm:

@jit
def DistanceEstimator(positions, iterations, degree=8, bailout=1000):
    m = positions.shape[0]
    x, y, z = np.zeros(m), np.zeros(m), np.zeros(m)
    x0, y0, z0 = positions[:, 0], positions[:, 1], positions[:, 2]
    dr = np.zeros(m) + 1
    r = np.zeros(m)
    theta = np.zeros(m)
    phi = np.zeros(m)
    zr = np.zeros(m)
    for _ in range(iterations):
        r = np.sqrt(x*x + y*y + z*z)
        idx1 = r < bailout
        dr[idx1] = np.power(r[idx1], degree - 1) * degree * dr[idx1] + 1.0

        theta[idx1] = np.arctan2(np.sqrt(x[idx1]*x[idx1] + y[idx1]*y[idx1]), z[idx1])
        phi[idx1] = np.arctan2(y[idx1], x[idx1])

        zr[idx1] = r[idx1] ** degree
        theta[idx1] = theta[idx1] * degree
        phi[idx1] = phi[idx1] * degree

        x[idx1] = zr[idx1] * np.sin(theta[idx1]) * np.cos(phi[idx1]) + x0[idx1]
        y[idx1] = zr[idx1] * np.sin(theta[idx1]) * np.sin(phi[idx1]) + y0[idx1]
        z[idx1] = zr[idx1] * np.cos(theta[idx1]) + z0[idx1]

    return 0.5 * np.log(r) * r / dr

Plot a Mandelbulb

Let's plot a classic Mandelbulb with $p=8$ and zoom in on the resulting image.

7th
8th

Shapes of Mandelbulb

We can take other values for $p$ and see how the resulting pictures differ.

9th
10th

You can play with zooming and constructing fractal sets in this Datalore workbook.

Header image source: A ray-traced image of the 3D Mandelbulb for the iteration v ↦ v8 + c, CC BY-SA 3.0

Show Comments