Fractal Generator

Daily procedural fractals — Python + NumPy

[← Back to Portfolio]

Overview

Every day, the site generates a new fractal image automatically. The current date is MD5-hashed to produce a deterministic seed, and that seed controls every decision: which of 28 fractal types to render, which pre-catalogued location to zoom into, and which of 20 color palettes to apply. The image is cached as a PNG so generation only happens once per day regardless of how many people load the page.

The 28 types span several fundamentally different mathematical families: escape-time fractals (Mandelbrot, Julia, Burning Ship, and many variants), convergence/Newton fractals, IFS/chaos-game fractals (Barnsley Fern, Sierpinski, Dragon Curve), and transcendental fractals using complex trig and exponentials. Each algorithm produces completely different visual geometry from the same NumPy vectorisation pattern.

Python NumPy Pillow Quart APScheduler

Code

// Seeding from the date
The date string is MD5-hashed and the first 8 hex digits are read as an integer. That single number then drives every subsequent random decision through np.random.seed(seed) — fractal type selection, zoom location, and palette choice are all a pure function of the date.
def date_to_seed(date_str):
    # e.g. "2026-03-30" -> deterministic integer unique to that day
    return int(hashlib.md5(date_str.encode()).hexdigest()[:8], 16)
// Fractal type dispatch — 28 types
The seed modulo the length of the type list picks a fractal family for the day. The same seed is passed to each generator, which uses it to select a zoom location and then calls np.random.seed(seed) for reproducible random variation within that location.
fractal_types = [
    'mandelbrot', 'julia', 'burning_ship', 'tricorn', 'multibrot',
    'newton', 'phoenix', 'celtic', 'cosine', 'exponential',
    'buffalo', 'perpendicular_burning_ship', 'mandelbar_julia',
    'nova', 'newton_z4', 'sine', 'manowar', 'glynn',
    'lambda', 'magnet1', 'barnsley_fern', 'sierpinski',
    'dragon_curve', 'collatz', 'hyperbolic_sine',
    'zubieta', 'burning_julia', 'mandelbrot_power',
]

def generate_fractal(date_str):
    seed = date_to_seed(date_str)
    fractal_type = fractal_types[seed % len(fractal_types)]

    # dispatch to the appropriate generator
    generators = {
        'mandelbrot': (generate_mandelbrot, 'Mandelbrot Set'),
        'julia':      (generate_julia,      'Julia Set'),
        'newton':     (generate_newton,     'Newton Fractal'),
        'phoenix':    (generate_phoenix,    'Phoenix Fractal'),
        'barnsley_fern': (generate_barnsley_fern, 'Barnsley Fern'),
        # ... 23 more
    }
    fn, name = generators[fractal_type]
    M = fn(seed)

    palette = generate_color_palette(seed)
    img_data = palette[M]                         # fancy index: each pixel gets its RGB
    img = Image.fromarray(img_data.astype('uint8'), 'RGB')
    img.save(os.path.join(CACHE_DIR, f"{date_str}.png"), 'PNG')
// Escape-time algorithm — Mandelbrot
The canonical approach: for each pixel c in the complex plane, iterate Z = Z² + c and record how many steps before |Z| > 2 (the escape radius). The iteration count becomes the pixel colour. Vectorised with NumPy so all 640,000 pixels run in parallel — the boolean mask stops updating pixels that have already escaped.
def generate_mandelbrot(seed):
    np.random.seed(seed)
    interesting_points = [
        (-0.7,     0.0,      1.5),    # main body
        (-0.5,     0.0,      0.5),    # seahorse valley
        (-0.16,    1.0405,   0.026),  # spiral
        (-0.7269,  0.1889,   0.005),  # double spiral
        (0.285,    0.01,     0.013),  # elephant valley
        (-0.748,   0.1,      0.014),  # triple spiral
        (-0.235125,0.827215, 0.004),  # satellite
        (-0.8,     0.156,    0.04),   # scepter valley
        # ... 4 more
    ]
    cx, cy, base_zoom = interesting_points[seed % len(interesting_points)]
    # small seeded nudge so the same location looks different on repeat visits
    cx += (np.random.rand() - 0.5) * base_zoom * 0.3
    cy += (np.random.rand() - 0.5) * base_zoom * 0.3
    zoom = base_zoom * (0.7 + np.random.rand() * 0.6)

    x = np.linspace(cx - zoom, cx + zoom, WIDTH)   # WIDTH = HEIGHT = 800
    y = np.linspace(cy - zoom, cy + zoom, HEIGHT)
    X, Y = np.meshgrid(x, y)
    C = X + 1j * Y                                # complex plane as a 2D array

    Z = np.zeros_like(C)
    M = np.zeros(C.shape, dtype=int)               # escape-time map

    for i in range(MAX_ITER):                       # MAX_ITER = 256
        mask = np.abs(Z) <= 2                      # pixels still inside escape radius
        Z[mask] = Z[mask]**2 + C[mask]             # Mandelbrot recurrence
        M[mask] = i                                 # overwrite until escaped

    return M
// Escape-time variants — changing one line of math
Most of the 28 types are variations on the same vectorised loop. The only thing that changes is the recurrence applied to Z on each step. The Burning Ship takes the absolute values of the real and imaginary parts before squaring, creating a jagged "ship" silhouette. The Celtic variant takes the absolute value of only the real part of Z², warping the set into a fractal knot pattern. The Phoenix fractal adds a feedback term using the previous iteration's Z value.
# Standard Mandelbrot:
Z[mask] = Z[mask]**2 + C[mask]

# Burning Ship — abs on both components before squaring:
Zr, Zi = Z[mask].real, Z[mask].imag
Z[mask] = (np.abs(Zr) + 1j * np.abs(Zi))**2 + C[mask]

# Celtic — abs on real part of Z² only:
Zr, Zi = Z[mask].real, Z[mask].imag
Z[mask] = (np.abs(Zr**2 - Zi**2) + 2j * Zr * Zi) + C[mask]

# Tricorn (Mandelbar) — conjugate of Z before squaring:
Z[mask] = np.conj(Z[mask])**2 + C[mask]

# Multibrot — higher powers (p = 3, 4, or 5):
Z[mask] = Z[mask]**p + C[mask]

# Phoenix — carries the previous Z as a feedback term:
Z_new = np.where(mask, Z**2 + c + p * Z_prev, Z)
Z_prev = np.where(mask, Z, Z_prev)
Z = Z_new

# Man o' War — C itself evolves alongside Z:
Z_new  = np.where(mask, Z**2 + C_cur, Z)
C_new  = np.where(mask, Z  + C_cur, C_cur)
Z, C_cur = Z_new, C_new
// Newton fractal — convergence, not escape
Newton fractals use a completely different criterion. Instead of checking whether Z escapes to infinity, each pixel applies Newton's root-finding method to z³ − 1 = 0, then checks which of the three cube roots of unity it converges to. The root index sets the colour region; the iteration count controls shading within that region. The boundary between basins of attraction is where the fractal structure lives.
def generate_newton(seed):
    np.random.seed(seed)
    zoom = 1.5 * (0.8 + np.random.rand() * 0.4)
    cx = (np.random.rand() - 0.5) * 0.4
    cy = (np.random.rand() - 0.5) * 0.4

    x = np.linspace(cx - zoom, cx + zoom, WIDTH)
    y = np.linspace(cy - zoom, cy + zoom, HEIGHT)
    X, Y = np.meshgrid(x, y)
    Z = X + 1j * Y

    # the three cube roots of unity: 1, e^(2πi/3), e^(4πi/3)
    roots = np.array([1, np.exp(2j * np.pi / 3), np.exp(4j * np.pi / 3)])
    root_map = np.full(Z.shape, -1, dtype=int)   # which root each pixel converges to
    iter_map = np.zeros(Z.shape, dtype=int)        # how many steps to convergence

    for i in range(MAX_ITER):
        denom = 3 * Z**2
        safe = np.abs(denom) > 1e-10              # avoid division by zero at critical pts
        # Newton step: z - f(z)/f'(z) = z - (z³-1)/(3z²) = (2z³+1)/(3z²)
        Z[safe] = Z[safe] - (Z[safe]**3 - 1) / denom[safe]

        for r_idx, root in enumerate(roots):
            converged = (np.abs(Z - root) < 1e-6) & (root_map == -1)
            root_map[converged] = r_idx
            iter_map[converged] = i

    # encode root basin + convergence speed into a single value for colouring
    M = root_map * (MAX_ITER // 3) + (iter_map % (MAX_ITER // 3))
    return np.clip(M, 0, MAX_ITER - 1)
// IFS / chaos game — Barnsley Fern
IFS (Iterated Function Systems) fractals work nothing like escape-time. Instead of testing every pixel, you run a single point through 500,000 random affine transformations drawn with specific probabilities, plotting each landing position. The Barnsley Fern uses four transforms: stem (1%), main frond (85%), left leaflet (7%), right leaflet (7%). The point cloud is then mapped to a pixel grid with log-scaling to handle density variation.
def generate_barnsley_fern(seed):
    np.random.seed(seed)
    tilt = (seed % 7) * 0.003    # subtle seed-based tilt variation

    # four affine transforms: matrix A, translation b, probability p
    transforms = [
        (np.array([[0, 0],    [0, 0.16]]),                        np.array([0, 0]),    0.01),  # stem
        (np.array([[0.85, 0.04+tilt], [-0.04+tilt, 0.85]]),     np.array([0, 1.6]),  0.85),  # main frond
        (np.array([[0.2, -0.26], [0.23, 0.22]]),               np.array([0, 1.6]),  0.07),  # left
        (np.array([[-0.15, 0.28],[0.26, 0.24]]),               np.array([0, 0.44]), 0.07),  # right
    ]
    probs = np.cumsum([t[2] for t in transforms])

    x, y = 0.0, 0.0
    n = 500_000
    xs, ys = np.empty(n), np.empty(n)

    for k in range(n):
        i = int(np.searchsorted(probs, np.random.rand()))
        A, b, _ = transforms[i]
        x, y = A @ np.array([x, y]) + b   # apply affine transform
        xs[k], ys[k] = x, y

    # map the point cloud onto a log-scaled density grid
    return _ifs_to_grid(xs, ys, flip_y=True)


def _ifs_to_grid(xs, ys, flip_y=True):
    # normalise coordinates to pixel indices
    xi = ((xs - xs.min()) / (xs.ptp() + 1e-9) * (WIDTH  - 1)).astype(int)
    yi = ((ys - ys.min()) / (ys.ptp() + 1e-9) * (HEIGHT - 1)).astype(int)
    if flip_y: yi = (HEIGHT - 1) - yi

    M = np.zeros((HEIGHT, WIDTH), dtype=float)
    np.add.at(M, (yi, xi), 1)              # accumulate hit counts

    M_log = np.log1p(M)                     # log-scale dampens density hotspots
    return (M_log / M_log.max() * (MAX_ITER - 1)).astype(int)
// Color palette — construction and application
Twenty named palettes (Fire, Ocean, Midnight, Matrix, etc.) are defined as 5–6 RGB anchor colours. The seed picks a palette, then NumPy linearly interpolates between the anchors to produce a 256-entry lookup table. Applying the palette to the iteration-count array is a single fancy-index operation — one line maps 640,000 integers to 640,000 RGB triples simultaneously.
def generate_color_palette(seed):
    np.random.seed(seed)
    palette_type = seed % 20

    # each entry is a list of RGB anchor colours
    palettes = {
        0:  # Fire
            [[0,0,0], [128,0,128], [255,0,0], [255,128,0], [255,255,0]],
        1:  # Ocean
            [[0,0,64], [0,64,128], [0,128,255], [128,200,255], [255,255,255]],
        12: # Matrix
            [[0,0,0], [0,32,0], [0,128,0], [0,255,0], [200,255,200]],
        # ... 17 more
    }
    colors = np.array(palettes.get(palette_type, palettes[0]))

    # interpolate anchors into a 256-entry lookup table
    palette = np.zeros((MAX_ITER, 3), dtype=np.uint8)
    for i in range(MAX_ITER):
        idx = (i / MAX_ITER) * (len(colors) - 1)
        lo, hi = int(idx), min(int(idx) + 1, len(colors) - 1)
        t = idx - lo
        palette[i] = colors[lo] * (1 - t) + colors[hi] * t

    return palette   # shape (256, 3)


# applying the palette: the entire 800x800 image in one line
img_data = palette[M]                              # (800,800) int -> (800,800,3) uint8
img = Image.fromarray(img_data.astype('uint8'), 'RGB')
// Transcendental fractals — complex trig and exp
Applying standard mathematical functions to Z in the complex plane produces completely different geometry from polynomial escape-time sets. cos(z) is entire (no poles), so the escape radius rises to 10 rather than 2. e^z maps vertical strips periodically, producing infinite repeating structures along the imaginary axis. The Collatz fractal uses a smooth analytic extension of the integer Collatz conjecture — the same recurrence that's never been proven to converge for all integers, visualised over the complex plane.
# Cosine fractal: Z = cos(Z) + C — escape radius 10 (trig grows faster than poly)
for i in range(MAX_ITER):
    mask = np.abs(Z) <= 10
    Z[mask] = np.cos(Z[mask]) + C[mask]
    M[mask] = i

# Exponential fractal: Z = e^Z + C — periodic along imaginary axis
for i in range(MAX_ITER):
    mask = np.abs(Z) <= 10
    Z[mask] = np.exp(Z[mask]) + C[mask]
    M[mask] = i

# Smooth Collatz extension over C: f(z) = (2 + 7z - (2 + 5z)·cos(πz)) / 4
# This is Hailstone iteration, analytically continued — identical to integer Collatz on ℤ
for i in range(MAX_ITER):
    mask = np.abs(Z) <= 10
    Z[mask] = (2 + 7*Z[mask] - (2 + 5*Z[mask]) * np.cos(np.pi * Z[mask])) / 4
    M[mask] = i
// Dragon Curve IFS — two-transform chaos game
The Heighway Dragon Curve is an IFS with just two transforms, each chosen with equal probability. Both are rotations-with-scale: the first rotates 45° around the origin, the second rotates 135° around the point (1, 0). Running 500,000 iterations of the chaos game traces out the self-similar dragon shape that would also emerge by folding a strip of paper in half repeatedly. The result goes through _ifs_to_grid — the same log-scaled density accumulator used for all IFS types.
# Two affine transforms for the Heighway dragon IFS
transforms = [
    (np.array([[ 0.5, -0.5], [ 0.5,  0.5]]), np.array([0, 0]), 0.5),  # rotate +45° around origin
    (np.array([[-0.5, -0.5], [ 0.5, -0.5]]), np.array([1, 0]), 0.5),  # rotate +135° around (1,0)
]
probs = np.cumsum([t[2] for t in transforms])

x, y = 0.0, 0.0
for k in range(500_000):
    i = int(np.searchsorted(probs, np.random.rand()))
    A, b, _ = transforms[i]
    x, y = A @ np.array([x, y]) + b    # apply chosen affine transform
    xs[k], ys[k] = x, y

return _ifs_to_grid(xs, ys, flip_y=False)
// Caching and serving
Generated images are stored as PNGs named by date (2026-03-30.png). The API checks the cache before generating, so repeat requests are just a file read. Metadata (fractal type, palette name, seed) is kept in a JSON sidecar. A Quart route serves the raw PNG bytes directly from disk.
def get_fractal(date_str):
    'Get cached fractal or generate it if this is the first request today.'
    filepath = os.path.join(CACHE_DIR, f"{date_str}.png")
    if os.path.exists(filepath):
        return filepath, load_metadata().get(date_str, {})
    return generate_fractal(date_str)     # generates, saves PNG + metadata, returns path


@fractal_bp.route('/api/today')
async def get_today():
    today = datetime.now().strftime('%Y-%m-%d')
    filepath, meta = get_fractal(today)
    return jsonify({
        'date':      today,
        'image_url': f'/fractal/api/image/{today}',
        'metadata':  meta        # {type, name, seed, palette, filename}
    })


@fractal_bp.route('/api/image/<date_str>')
async def get_image(date_str):
    filepath, _ = get_fractal(date_str)
    return await send_file(filepath, mimetype='image/png')

Links

[View Today's Fractal] [Technical Paper →] [← Portfolio]