Mandelbox fold-scale 2D projection

Math essence

For each c, iterate z := boxFold(z); if |z|<0.5 then z := 4z else if |z|<1 then z := z/|z|^2; then z := scale*z + c. Escape when |z| > 2.5. boxFold mirrors each component z_i -> 2*clamp(z_i, -1, 1) - z_i, producing a [0, 1] -> [-1, 0] reflection. The (u, v) grid is rotated by pitch/yaw and averaged over 2 z-planes, then coloured by minimum |z| along the orbit and by the iteration count of escape.

Parameters

Reference

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

from __future__ import annotations

from typing import ClassVar

import numpy as np


def _mandelbox_iterate(
    cx: np.ndarray, cy: np.ndarray, cz: np.ndarray, max_iter: int, escape: float,
    scale: float, fold: float
) -> tuple[np.ndarray, np.ndarray]:
    """Tom Lowe's Mandelbox iteration (2010).

    For each starting point c = (cx, cy, cz), iterate
        z := boxFold(z)
        if |z| < 0.5:   z := z * 4
        elif |z| < 1:   z := z / |z|^2
        z := scale * z + c
    Escape when |z| > escape.

    The boxFold is the per-component fold used in IQ/MIQEL's reference
    implementations:  z_i := clamp(z_i, -1, 1) * 2 - z_i, which mirrors
    the cube [-1, 1]^3 onto itself.  We track the minimum |z| (orbit
    trap) along each trajectory, the same quantity that gives the
    dramatic colour bands in scale-2 and scale-3 Mandelbox renders.
    """
    x = cx.copy()
    y = cy.copy()
    z = cz.copy()
    iters = np.zeros(cx.shape, dtype=np.int32)
    orbit_min = 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
        r_m = np.sqrt(r2[mask])
        orbit_min[mask] = np.minimum(orbit_min[mask], r_m.astype(np.float32))

        # Box-fold: mirror each component across the [-1, 1] range
        xs = np.clip(x[mask], -fold, fold)
        ys = np.clip(y[mask], -fold, fold)
        zs = np.clip(z[mask], -fold, fold)
        xs = xs * 2.0 - x[mask]
        ys = ys * 2.0 - y[mask]
        zs = zs * 2.0 - z[mask]

        # Sphere-fold
        r2_s = xs * xs + ys * ys + zs * zs
        r_s = np.sqrt(r2_s)
        small = r_s < 0.5
        med = (r_s >= 0.5) & (r_s < 1.0)
        if np.any(small):
            xs[small] *= 4.0
            ys[small] *= 4.0
            zs[small] *= 4.0
        if np.any(med):
            r2_m = r2_s[med]