The Icebox Fractal

The Mandelbulb and Mandelbox are fascinating 3D fractals, but they lack a satisfying representation of my favorite fractal: ice. They have swirls, bulbs, mysterious tunnels, and half-constructed skyscrapers, but none of the elaborate layered crystals one sees in a single snowflake.

So, I created a new fractal specifically to resemble ice, by combining and modifying parts of those two fractal equations. The Mandelbulb has a sweeping natural surface with excellent symmetry, but the smooth rotations in the equation cannot have sharp geometry. I introduced sharp fractured features by modifying the polar coordinates to use a non-Euclidean rectangular distance metric that has a discontinuous first derivative. Instead of using x² + y² + z² as the radius, I used max(|x|, |y|, |z|), which picks the largest distance projected along one of the basis vectors. Geometrically, this distance function represents a square pyramid instead of a cone.

The sharp creases of the pyramid are where the first derivative is not continuous. The abrupt change in direction introduces sharp features into the fractal when the iterated point crosses those boundaries. The rectangular shape of the pyramid makes the new fractal shape a box instead of a bulb, since in this distance metric, points on the surface of a cube are all the same distance from the origin. The result is a quite satisfying icy box!

However, there is one final improvement. The Mandelbox's folding calculations are relatively generic for fractals, and can be applied as an additional effect. After the polar coordinate rotation, the rectangular coordinates are then cube-folded and sphere-inverted using parts of the Mandelbox equation. The inversion hollows out the inside, and dramatically enhances the frosty features!

(Click to open in the Icebox viewer.)

The cube folding creates a solid effect. At small values, it reveals icy swirls from the rotations in the equation.

(Click to open in the Icebox viewer.)

At large values, the folding fractures the box into infinite nested ice cubes.

(Click to open in the Icebox viewer.)

Finally, as with most fractals, the details become exponentially more intracate as you zoom in!

(Click to open in the Icebox viewer.) Ice Crystals in the Icebox Fractal

(Click to open in the Icebox viewer.)

All these beautiful ice images, and infinitely many more, are contained in this simple equation. (This code is presented as a WebGL distance estimator function, which is how the viewer renders the fractal in real time on the GPU. It returns the distance from a test point to the fractal boundary. The fractal set equation is slightly simpler. It skips the calculation of dr, and just returns false if r exceeds 2, otherwise true.)

highp float icebox_distance_estimator(highp vec3 c, highp float mPower, highp float mFoldLimit) {
    highp vec3 z = c;
    highp float r, r2;
    highp float dr = 1.0;
    
    for (int n = 0; n < 32; n++) { // Use more iterations to resolve fine features at high zoom levels.
        r = max(abs(z.x), max(abs(z.y), abs(z.z)));
        if (r > 2.0) break;
        
        // Use polar coordinates. (Calculate the Mandelbulb equation.)
        highp float theta = acos(z.z/r);
        highp float phi = atan(z.y, z.x);
        highp float derivpow = pow(r, mPower - 1.0);
        dr = derivpow*mPower*dr + 1.0;
        
        // Calculate the generalized Mandelbrot equation: z = z^k + c
        theta = theta*mPower; // Complex exponentiation represents multiplication of the vector angle.
        phi = phi*mPower;
        
        // Convert back to rectangular coordinates, so we can add the sample point.
        z = (derivpow*r)*vec3(cos(phi)*sin(theta), sin(phi)*sin(theta), cos(theta));
        z += c;
        
        // Apply box folds and sphere scaling from the Mandelbox equation.
        z.xyz = clamp(z.xyz, -mFoldLimit, mFoldLimit)*2.0 - z.xyz;
        r2 = dot(z.xyz, z.xyz);
        if (r2 < 0.5) { // Inner Scaling
            z *= 2.6;
        } else if (r2 < 1.3) { // Sphere Inversion
            z *= (1.3/r2);
        }
    }
    return 0.25*log(r)*r/dr; // Return the distance estimate from c to the fractal boundary.
}

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.