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

#### 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.

#### 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
```

#### 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
```

## 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.

### Shapes of Mandelbulb

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

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*