Previous posts: part I, part II, and part III.

Despite its young age, the Mandelbulb is probably the most famous 3D fractal in existence. This post will examine how we can create a Distance Estimator for it. But before we get to the Mandelbulb, we will have to step back and review a bit of the history behind it.

The Search for the Holy Grail

The original Mandelbrot fractal is a two dimensional fractal based on the convergence properties of a series of complex numbers. The formula is very simple: for any complex number z, check whether the sequence iteratively defined by \(z_{n+1} = z_{n}^2+c\) diverges or not. The Mandelbrot set is defined as the set of points which do not diverge, that is, the points with a series that stays bounded within a given radius. The results can be depicted in the complex plane.

The question is how to extend this to three dimensions. The Mandelbrot set fits two dimensions, because complex numbers have two components. Can we find a similar number system for three dimensions?

The Mandelbrot formula involves two operations: adding numbers, and squaring them. Creating a n-component number where addition is possible is easy. This is what mathematicians refer to as a vector space. Component-wise addition will do the trick, and seems like the logical choice.

But the Mandelbrot formula also involves squaring a number, which requires a multiplication operator (a vector product) on the vector space. A vector space with a (bilinear) vector product is called an algebra over a field. The numbers in these kind of vector spaces are often called hypercomplex numbers.

To see why a three dimensional number system might be problematic, let us try creating one. We could do this by starting out with the complex numbers and introduce a third component, j. We will try to keep as many as possible of the characteristic properties of the complex and real numbers, such as distributivity, \(a*(b+c)=(a*b)+(a*c)\), and commutativity, \(a*b=b*a\). If we assume distributivity, we only need to specify how the units of the three components multiply. This can be illustrated in a multiplication table. Since we also assumed commutativity, such a table must be symmetric:

\(
\begin{bmatrix}
& \boldsymbol{1} & \boldsymbol{i} & \boldsymbol{j} \\
\boldsymbol{1} & 1 & i & j \\
\boldsymbol{i} & i & -1 & ?\\
\boldsymbol{j} & j & ? & ?
\end{bmatrix}
\)
 
For a well-behaved number system, anything multiplied by 1 is still one, and if we now require the real and imaginary components to behave as for the complex numbers, we only have three components left – the question marks in the matrix. I’ve rendered out a few of the systems, I encountered while trying arbitrary choices of the missing numbers in the matrix:

(Many people have explored various 3D component multiplication tables – see for instance Paul Nylander’s Hypercomplex systems for more examples).

Unfortunately, our toy system above fails to be associative (i.e. it is not always true, that \(a*(b*c) = (a*b)*c\)), as can be seen by looking at the equation \(i*(i*j) = (i*i)*j \Rightarrow i*x = -j\), which can not be satisfied no matter how we choose x.

It turns out that it is difficult to create a consistent number system in three dimensions. There simply is no natural choice. In fact, if we required that our number system allowed for a division operator, there is a theorem stating that only four such mathematical spaces are possible: the real numbers (1D), the complex numbers (2D), the quaternions (4D) and the octonions (8D). But no 3D systems.

But what about the 4D Quaternions? Back in 1982, Alan Norton published a paper showing a Quaternion Julia set made by displaying a 3D “slice” of the 4D space. Here is an example of a Quaternion Julia fractal:

Of course, in order to visualize a 4D object, you have to make some kind of dimensional reduction. The most common approach is to make a 3D cross-section, by simply keeping one of the four components at a fixed value.

If you wonder why you never see a Quaternion Mandelbrot image, the reason is simple. It is not very interesting because of its axial symmetry:

If you, however, make a rotation inside the iteration loop, you can get something more like a 3D Mandelbrot.

The Quaternion system (and the 3D hypercomplex systems above) are defined exactly as the 2D system – by checking if \(z_{n+1} = z_{n}^2+c\) converges or not.

But how do we draw a 3D image of these fractals? In contrast to the 2D case, where it is possible to build a 2D grid, and check inside each cell, building a 3D grid and checking each cell would be far too memory and time consuming for images in any decent resolution.

A distance estimator for quadratic systems.

While Alan Norton used a different rendering approach, a very elegant solution to this was found by John Hart et al in a 1989 paper: distance estimated rendering. As discussed in the previous posts, distance estimated rendering requires that we are able to calculate a lower bound to the distance from every point in space to our fractal surface! A first, this might seem impossible. But it turns out such a formula already was known for the 2D Mandelbrot set. A distance estimate can be found as:

(1)   \(DE=0.5*ln(r)*r/dr\)
 
Where ‘r’ is the escape time length, and ‘dr’ is the length of the running derivative. (The approximation is only exact in the limit where the number of iterations goes to infinity)

In order to define what we mean by the running derivative, we need a few extra definitions. For Mandelbrot sets, we study the sequence \(z_{n+1} = z_{n}^2+c\) for each point c. Let us introduce the function \(f_n(c)\), defined as the n’th entry for the sequence for the point c. By this definition, we have the following defining formula for the Mandelbrot set:

\(f_n(c) = f_{n-1}^2(c) + c, f_0(c) = 0\)
 
Deriving this function with respect to c, gives

(2)   \(f’_n(c) = 2f_{n-1}(c)f’_{n-1}(c)+1\)   (for Mandelbrot formula)
 
Similar, the Julia set is defined by choosing a fixed constant, d, in the quadratic formula, using c only as the first entry in our sequence:

\(f_n(c) = f_{n-1}^2(c) + d, f_0(c) = c\)
 
Deriving this function with respect to c, gives

(3)   \(f’_n(c) = 2f_{n-1}(c)f’_{n-1}(c)\)   (for Julia set formula)
 
which is almost the same result as for the Mandelbrot set, except for unit term. And now we can define the length of \(f_n\), and the running derivative \(f’_n\):
 
\(r = |f_n(c)|\) and \(dr = |f’_n(c)|\)
 
used in the formula (1) above. This formula was found by Douady and Hubbard in a 1982 paper (more info).

2D Julia set rendered using a distance estimator approach. This makes it possible to emphasize details, without having to use extensive oversampling.

Due to a constraint in WordPress, this post has reached its maximum length. The next post continues the discussion, and shows how the formula above can be used for other types of fractals than the 2D Mandelbrot.