Mandelbulb power-8 2D projection

Math essence

In spherical coordinates (r, theta, phi) Quilez's iteration is r' = r^8, theta' = 8 theta, phi' = 8 phi followed by w' = w' + c. The 2D projection samples the field on a regular (u, v) grid, applies a pitch/yaw rotation and zoom, averages 3 z-planes, and mixes the per-axis orbit-trap minima with a slow colour ramp so the surface tessellation of the bulb is visible. step(n) morphs camera_pitch/yaw/zoom sinusoidally.

Parameters

Reference

First 60 lines of reference.py. Full source: library_external/fractals_3d_webgl/Mandelbulb2DProjection/reference.py

from __future__ import annotations

from typing import ClassVar

import numpy as np


def _mandelbulb_pow8_iterate(
    cx: np.ndarray, cy: np.ndarray, cz: np.ndarray, max_iter: int, escape: float
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Quilez's power-8 Mandelbulb iteration, fully vectorized.

    For each starting point c = (cx, cy, cz), iterate w = w^8 + c using
    Quilez's equal-length-exponentiation / equal-angle-multiplication scheme
    (iquilezles.org/articles/mandelbulb). Concretely: in spherical
    coordinates (wr, wo=acos(y/r), wi=atan2(x,z)) update
        wr <- wr^8
        wo <- 8 * wo
        wi <- 8 * wi
    and convert back to Cartesian, then add c. This is *not* a quaternionic
    square (it is intentionally algebraically inconsistent) but produces the
    iconic bulb geometry.

    Returns (iters, ox, oy, oz): iteration count to escape, plus the
    *minimum* (orbit-trap) values of |w_x|, |w_y|, |w_z| seen during the
    orbit. Those three orbit traps are the same ones Quilez uses to colour
    the bulb: closer to a coordinate plane => lower trap value => a
    different "region" of the surface.
    """
    x = cx.copy()
    y = cy.copy()
    z = cz.copy()
    iters = np.zeros(cx.shape, dtype=np.int32)
    ox = np.full(cx.shape, np.inf, dtype=np.float32)
    oy = np.full(cx.shape, np.inf, dtype=np.float32)
    oz = np.full(cx.shape, np.inf, dtype=np.float32)
    escape2 = escape * escape
    for _ in range(max_iter):
        r2 = x * x + y * y + z * z
        mask = r2 < escape2
        if not np.any(mask):
            break
        iters[mask] += 1
        ox = np.minimum(ox, np.abs(x).astype(np.float32))
        oy = np.minimum(oy, np.abs(y).astype(np.float32))
        oz = np.minimum(oz, np.abs(z).astype(np.float32))
        r = np.sqrt(r2[mask])
        # Spherical update
        wr = r ** 8
        theta = np.arccos(np.clip(y[mask] / r, -1.0, 1.0)) * 8.0
        phi = np.arctan2(x[mask], z[mask]) * 8.0
        sin_t = np.sin(theta)
        nx = wr * sin_t * np.sin(phi) + cx[mask]
        ny = wr * np.cos(theta) + cy[mask]
        nz = wr * sin_t * np.cos(phi) + cz[mask]
        x[mask] = nx
        y[mask] = ny
        z[mask] = nz
    ox = np.where(np.isfinite(ox), ox, 1.0)
    oy = np.where(np.isfinite(oy), oy, 1.0)