The Mandelbulb
Fractals can exist in any number of dimensions. The classic Mandelbrot set can be generalized to 3D and beyond. In its original 2D form, the Mandelbrot set is defined as the set of all complex numbers, C, where if you repeat the following formula forever, Z does not grow infinitely large:
Z = Z2 + C
But what makes this special? Geometrically, it is a squaring and a rotation. A complex number is just a point in two dimensions. It can be expressed in polar coordinates as a direction, θ, and distance, r, outward from the origin.
Using Euler's Formula, the distance and direction can be written directly as a complex number in exponential notation:
z = reiθ = r(cos(θ) + isin(θ))Therefore, using elementary algebra:
z2 = (reiθ)(reiθ) = r2eiθ + iθ = r2e2iθIn other words, squaring a complex number squares the distance and rotates the point by its own angle.
These rotations are what create the infinite swirly patterns in the Mandelbrot set. For each sample point, it rotates around and around by repeatedly adding the point, C, and squaring the cumulative value, Z, until it either diverges to infinity or stays finite. However, the exact angle of rotation depends on the point itself, which creates different harmonics for even the most subtle changes, that add up over many iterations to produce wildly different patterns!

↗ Click to open in the Mandelbrot viewer.
But how do we extend this to three dimensions? Complex numbers have no equivalent in three dimensions. (The "next" set of numbers is quaternions, which exist in 4D. The Mandelbrot formula can be calculated in 4D and projected into 3D. Unfortunately, the fractal is mostly a blob of whipped cream, and not the majestic maze of swirls and curls it is in 2D.)
Instead, we use spherical coordinates. A point in 3D can be represented as a latitude and longitude on a sphere (like the Earth), where the radius of the sphere is the distance from the origin*. Instead of (x, y, z), the three components are (ρ, φ, θ), where those letters (rho, phi, theta) represent (radius, latitude, longitude).
To extend the Mandelbrot formula to 3D, the cumulative value, Z, is converted to spherical coordinates in each iteration, the radius is squared, and both angles are doubled: (ρ2, 2φ, 2θ) Then it is converted back to rectangular coordinates and the sample point, C, is added. On the xy-plane, the result is identical to the Mandelbrot set, since the zero latitude remains unchanged after multiplication. But there are now extra swirls that extend outward in three dimensions, which you can see if you rotate the shape in the viewer.

↗ Click to open in the Mandelbulb viewer.
However, this is more of a Mandelbrot pancake than a 3D fractal! To increase the interesting details, we can apply a constant offset to the polar angle φ to extend the swirls further into 3D space:

↗ Click to open in the Mandelbulb viewer.
The rotation is the key part of the generalized formula that produces Mandelbrot-like swirls. The power is actually an arbitrary choice in the Mandelbrot set equation:
Zn = Znk + C = rkeikθ + C for any k ∈ ℤ ≥ 2In the formula for 3D, there is a radius and two angle components, but the same procedure applies.
The Mandelbulb Set Equation
Exponentiate the radius by k and multiply the angles by k: (ρk, kφ, kθ) at each iteration, where the point to test C = (x, y, z), and Z0 = (0, 0, 0), for C, Zn ∈ ℝ3
Let ρ = |Zn|, φ = sin−1(Then Zn+1 = ρk(cos(kθ)cos(kφ), sin(kθ)cos(kφ), sin(kφ)) + C
If the magnitude of Zn diverges, then the point C is not inside the Mandelbulb. Otherwise, if it stays bounded, then the point is inside. Although the squared equation successfully preserves the shape of the Mandelbrot set, the extra dimension introduces a lot of wasted space between the swirls. This is where the arbitrary power k comes in. When you raise Z to a higher power, it creates more bulb-like features in the fractal, hence the name "Mandelbulb". Here is an image where k = 8, which is the most popular power.

↗ Click to open in the Mandelbulb viewer.
The polar angle φ can be animated with an offset to produce an effect where it continuously unfolds from within.
As with all fractals, the details truly shine when you zoom in!
↗ Click to open in the Mandelbulb viewer.
Other powers are interesting too, and you can experiment with them in the viewer.
It uses a distance estimator to render the fractal in 3D by marching rays toward it from the viewing position. It repeatedly calculates how far each ray is from the fractal set (using the distance estimator), and then advances the ray by that amount. It will converge on the surface of the fractal set eventually, and much faster than sampling points to see if they are inside or outside. Here is a GLSL function to calculate the distance from a point, c, to the surface of the Mandelbulb:
float distance_estimator(vec3 c, float polarOffset, float power) {
vec3 z = c;
float r;
float dr = 1.0;
for (int n = 0; n < 12; n++) { // Number of iterations.
r = length(z);
if (r > 2.0) break; // Stop if the point diverges.
// Use spherical coordinates.
float phi = asin(z.z/r) + polarOffset; // Offset to make an "unfolding" animated effect.
float theta = atan(z.y, z.x);
float derivpow = pow(r, power - 1.0);
dr = derivpow*power*dr + 1.0; // Estimate the derivative, for the Hubbard-Douady potential.
// Calculate the generalized Mandelbrot equation: z = z^k + c
phi = phi*power; // Complex exponentiation represents multiplication of the vector angle.
theta = theta*power;
// Convert back to rectangular coordinates and add the sample point c.
z = (derivpow*r)*vec3(cos(theta)*cos(phi), sin(theta)*cos(phi), sin(phi));
z += c;
}
return 0.5*log(r)*r/dr; // Return an estimate of the distance to the fractal set.
}
* I am using the geographic convention for polar coordinates, where the polar angle starts at 0° for points on the xy-plane, ranging between ± 90° as a measure of elevation (latitude). This allows the Mandelbrot set to remain unchanged on the xy-plane for a power-2 Mandelbulb, since the latitude angle of 0° is not affected by the angle multiplication. Other spherical coordinate conventions produce the same fractal, but require a constant offset for the polar angle to account for the differences.
All code snippets in this article are open source under the MIT license. You are free to use them for both personal and business uses.
