Penrose Diagram v3

I took another stab at generating a Penrose diagram using Chatgpt 5. The resulting Python script was executed in a Google Colab notebook:

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Penrose (conformal) diagram using every 24-bit RGB color exactly once.

This version includes:
  • STIX Two Math auto-download & use for all labels (covers ? and en dashes).
  • Superscripts (+, −, 0) composed with font metrics (no clipping, no tofu).
  • MANY semantic labels (?±, i±/i0, H± along horizons, ER bridge, etc.).
  • Crisp dark outline around the Penrose diamond (|x'| + |y'| = 1).
  • EVEN, matched r=const curves in all four regions (II/IV mirror I/III).
  • Center-bright radial gradients per region + low-chroma ripple background.
  • Sub-pixel horizon margin (no offsets).
  • Strict exclusive layer partitioning ⇒ exact all-RGB accounting.

Maximally extended Schwarzschild spacetime (Kruskal–Szekeres), units G=c=1, M=1.
"""

import os, time, math, pathlib
from datetime import datetime
import numpy as np
from PIL import Image, ImageDraw, ImageFont, ImageFilter, Image as PILImage

# Optional network fetch for font (works in Colab)
try:
    import requests
except Exception:
    requests = None  # graceful fallback if blocked

# ----------------------------
# Canvas & appearance
# ----------------------------
W = H = 4096
MARGIN_FRAC = 0.06
LINE_PX = 4             # horizons/singularities thickness
GRID_LINE_PX = 2
FRAME_PX = 6            # dark outline around the diamond
SING_MARGIN_FRAC = 0.08
LABEL_FONT_SIZE = 48
LABEL_STROKE_PIX = 2
VERIFY_UNIQUENESS = True
OUTPUT_DIR = "."

# ----------------------------
# Grid controls
# ----------------------------
DRAW_GRID = True
GRID_SAMPLES = 1700
KRUSKAL_K = 0.60

# r=const families
N_R_EXT = 14             # r>2 (Regions I & III)
N_R_INT = N_R_EXT        # r<2 (Regions II & IV), must match exterior
Y_REF_INT = 0.50         # interior reference level (|y'|), mid-wedge for robust spread

# Constant-t lines
N_T = 13
THETA_MAX_DEG = 42.0

# Optional non-radial null geodesics (OFF by default)
DRAW_CURVED_NULL = False
GEO_CURVE_PX = 2
GEODESIC_IMPACT_PARAMS = [5.6, 6.5, 8.0, 11.0, 16.0]
GEODESIC_RMAX = 60.0
GEODESIC_SAMPLES = 1800
COMPACT_K_NULL = 0.40

# Region hue personalities
REGION_HUE_PREFS = {
    "I":   [(210, 300)],          # blue–violet
    "III": [(20,  70)],           # orange–yellow
    "II":  [(300, 360), (0, 20)], # magenta–red
    "IV":  [(70,  200)],          # green–cyan
}
MIN_CHROMA_FOR_REGIONS = 25.0

# ----------------------------
# sRGB → LAB helpers (D65)
# ----------------------------
def srgb_to_linear(c):
    a = 0.055
    return np.where(c <= 0.04045, c/12.92, ((c + a) / (1 + a))**2.4)

def rgb01_to_lab(rgb01):
    M = np.array([[0.4124564, 0.3575761, 0.1804375],
                  [0.2126729, 0.7151522, 0.0721750],
                  [0.0193339, 0.1191920, 0.9503041]])
    Xn, Yn, Zn = 0.95047, 1.0, 1.08883
    lin = srgb_to_linear(rgb01)
    XYZ = np.dot(lin, M.T)
    X = XYZ[...,0] / Xn
    Y = XYZ[...,1] / Yn
    Z = XYZ[...,2] / Zn
    eps = (6/29)**3
    kappa = (29/3)**2/3
    def f(t): return np.where(t > eps, np.cbrt(t), (kappa*t + 16)/116)
    fX, fY, fZ = f(X), f(Y), f(Z)
    L = 116*fY - 16
    a = 500*(fX - fY)
    b = 200*(fY - fZ)
    return L, a, b

def atan2_deg(y, x):
    ang = np.degrees(np.arctan2(y, x))
    return np.where(ang < 0, ang + 360.0, ang)

# ----------------------------
# Palette allocator (64^3 bins; 64 colors/bin)
# ----------------------------
class BinPalette:
    def __init__(self):
        rb = np.arange(64, dtype=np.uint16)
        gb = np.arange(64, dtype=np.uint16)
        bb = np.arange(64, dtype=np.uint16)
        Rb, Gb, Bb = np.meshgrid(rb, gb, bb, indexing='ij')

        Rc = (Rb*4 + 2).astype(np.float32)
        Gc = (Gb*4 + 2).astype(np.float32)
        Bc = (Bb*4 + 2).astype(np.float32)
        rgb01 = np.stack([Rc, Gc, Bc], axis=-1) / 255.0

        L, a, b = rgb01_to_lab(rgb01)
        C = np.sqrt(a*a + b*b)
        h = atan2_deg(b, a)

        self.L = L.ravel(); self.C = C.ravel(); self.h = h.ravel()
        self.bin_count = 64**3
        self.consumed = np.zeros(self.bin_count, dtype=np.uint8)

        # intra-bin sequence (Morton-like) to spread subcolors
        offs = []
        for dr in range(4):
            for dg in range(4):
                for db in range(4):
                    morton = ((dr & 1) | ((dg & 1) << 1) | ((db & 1) << 2)
                              | (((dr>>1)&1) << 3) | (((dg>>1)&1) << 4) | (((db>>1)&1) << 5))
                    offs.append((morton, dr, dg, db))
        offs.sort(key=lambda t: t[0])
        self.intra = np.array([[dr,dg,db] for _,dr,dg,db in offs], dtype=np.uint8)

        self._rb = np.repeat(np.arange(64, dtype=np.uint16), 64*64)
        self._gb = np.tile(np.repeat(np.arange(64, dtype=np.uint16), 64), 64)
        self._bb = np.tile(np.arange(64, dtype=np.uint16), 64*64)

        self._dark_bins  = np.argsort(self.L)
        self._light_bins = np.argsort(-self.L)
        self._gray_bins  = np.lexsort((self.L, self.C))

    def _bin_base_rgb(self, bin_id):
        return np.array([int(self._rb[bin_id])*4,
                         int(self._gb[bin_id])*4,
                         int(self._bb[bin_id])*4], dtype=np.uint16)

    def _take_from_bin(self, bin_id, n):
        used = int(self.consumed[bin_id])
        if used >= 64: return None
        take = min(n, 64 - used)
        base = self._bin_base_rgb(bin_id)
        sl = self.intra[used:used+take].astype(np.uint16)
        cols = (base[None,:] + sl).astype(np.uint8)
        self.consumed[bin_id] = used + take
        return cols

    def allocate_from_bins(self, order, count):
        out = np.empty((count,3), dtype=np.uint8); filled = 0
        for bin_id in order:
            if filled >= count: break
            got = self._take_from_bin(bin_id, count - filled)
            if got is None: continue
            n = got.shape[0]; out[filled:filled+n] = got; filled += n
        if filled != count:
            for bin_id in range(self.bin_count):
                if filled >= count: break
                got = self._take_from_bin(bin_id, count - filled)
                if got is None: continue
                n = got.shape[0]; out[filled:filled+n] = got; filled += n
        assert filled == count, f"Allocator shortage: {filled}/{count}"
        return out

    def allocate_dark(self, count):       return self.allocate_from_bins(self._dark_bins,  count)
    def allocate_light(self, count):      return self.allocate_from_bins(self._light_bins, count)
    def allocate_low_chroma(self, count): return self.allocate_from_bins(self._gray_bins,  count)

    def bins_by_hue_preference(self, ranges, min_chroma=0.0):
        h, C, L = self.h, self.C, self.L
        def dist_to_ranges(hv):
            best = 180.0
            for a,b in ranges:
                if a<=b:
                    if a<=hv<=b: return 0.0
                    best = min(best, min(abs(hv-a), abs(hv-b)))
                else:
                    if (hv>=a) or (hv<=b): return 0.0
                    d1 = (a - hv) % 360.0; d2 = (hv - b) % 360.0
                    best = min(best, min(d1, d2))
            return best
        d = np.vectorize(dist_to_ranges)(h)
        order = np.lexsort((np.abs(L-50.0), -C, d))
        if min_chroma>0:
            hi = np.where(C[order] >= min_chroma)[0]
            lo = np.where(C[order] <  min_chroma)[0]
            order = np.concatenate([order[hi], order[lo]])
        return order

# ----------------------------
# Geometry & masks (diamond + regions + horizons)
# ----------------------------
def make_coordinate_grids(w, h, margin_frac):
    y, x = np.mgrid[0:h, 0:w].astype(np.float32)
    cx, cy = (w - 1) / 2.0, (h - 1) / 2.0
    half = min(w, h) * (0.5 - margin_frac)
    xp = (x - cx) / half
    yp = (cy - y) / half
    return xp, yp, half, cx, cy

def make_masks(w, h, margin_frac, line_px, sing_margin_frac):
    xp, yp, half, cx, cy = make_coordinate_grids(w, h, margin_frac)
    diamond = (np.abs(xp) + np.abs(yp) <= 1.0)
    s1 = yp - xp; s2 = yp + xp
    t = line_px / half

    horizons = ((np.abs(s1) <= t) | (np.abs(s2) <= t)) & diamond

    reg_II  = (s1 >  t) & (s2 >  t) & diamond
    reg_IV  = (s1 < -t) & (s2 < -t) & diamond
    reg_I   = (s1 < -t) & (s2 >  t) & diamond
    reg_III = (s1 >  t) & (s2 <  t) & diamond

    y_top =  1.0 - sing_margin_frac
    y_bot = -1.0 + sing_margin_frac
    singular = ((np.abs(yp - y_top) <= t) | (np.abs(yp - y_bot) <= t)) & diamond

    outside = ~diamond

    return {
        "diamond": diamond,
        "I": reg_I, "II": reg_II, "III": reg_III, "IV": reg_IV,
        "horizons": horizons, "singular": singular, "outside": outside,
        "xp": xp, "yp": yp, "half": half, "cx": cx, "cy": cy
    }

# ----------------------------
# Font: STIX Two Math (auto-download) + robust fallbacks
# ----------------------------
STIX_MATH_URLS = [
    # CTAN mirrors hosting the OTF
    "https://mirrors.ibiblio.org/pub/mirrors/CTAN/fonts/stix2-otf/STIXTwoMath-Regular.otf",
    "https://ftp.jaist.ac.jp/pub/CTAN/fonts/stix2-otf/STIXTwoMath-Regular.otf",
    "https://mirrors.ctan.org/fonts/stix2-otf/otf/STIXTwoMath-Regular.otf",
]

def ensure_stix_two_math_otf(local_dir="fonts", filename="STIXTwoMath-Regular.otf"):
    """Download STIX Two Math OTF if missing; return filesystem path."""
    pdir = pathlib.Path(local_dir); pdir.mkdir(parents=True, exist_ok=True)
    fpath = pdir / filename
    if fpath.exists() and fpath.stat().st_size > 0:
        return str(fpath)
    if requests is None:
        raise RuntimeError("requests not available to fetch STIXTwoMath-Regular.otf")
    last_err = None
    for url in STIX_MATH_URLS:
        try:
            r = requests.get(url, stream=True, timeout=30)
            r.raise_for_status()
            with open(fpath, "wb") as out:
                for chunk in r.iter_content(1 << 14):
                    if chunk: out.write(chunk)
            if fpath.stat().st_size > 0:
                return str(fpath)
        except Exception as e:
            last_err = e
    raise RuntimeError(f"Could not fetch STIXTwoMath-Regular.otf. Last error: {last_err}")

def load_font_main(size):
    """Prefer STIX Two Math (download if needed), then robust fallbacks."""
    try:
        path = ensure_stix_two_math_otf()
        return ImageFont.truetype(path, size)
    except Exception:
        for p in [
            "/usr/share/fonts/opentype/stix-two/STIXTwoMath-Regular.otf",
            "/usr/share/fonts/truetype/noto/NotoSansMath-Regular.ttf",
            "/usr/share/fonts/opentype/xits-math/XITSMath-Regular.otf",
            "/usr/share/fonts/opentype/lm-math/LatinModernMath-Regular.otf",
            "/usr/share/fonts/truetype/freefont/FreeSerif.ttf",
            "/usr/share/fonts/truetype/dejavu/DejaVuSerif.ttf",
            "DejaVuSerif.ttf",
        ]:
            try:
                return ImageFont.truetype(p, size)
            except Exception:
                continue
        return ImageFont.load_default()

# ----------------------------
# Labels → fill & stroke masks (rotations + composed superscripts)
# ----------------------------
def place_labels_masks(w, h, masks, font_size=48, stroke_px=2):
    """Return (fill_mask, stroke_mask) as boolean arrays, using STIX Two Math."""
    img = Image.new("L", (w, h), 0)
    drw = ImageDraw.Draw(img)
    font_main = load_font_main(font_size)
    cx, cy, half = masks["cx"], masks["cy"], masks["half"]

    def xy_from_norm(xp, yp):
        return (int(round(cx + xp * half)), int(round(cy - yp * half)))

    def draw_text(text, xp, yp, anchor="mm", angle=0, font=None):
        """Plain run (no superscript composition)."""
        if font is None: font = font_main
        if angle == 0:
            drw.text(xy_from_norm(xp, yp), text, fill=255, font=font, anchor=anchor)
        else:
            pad = max(12, int(0.45 * font.size))
            tw = drw.textlength(text, font=font)
            a, d = font.getmetrics()
            th = a + d
            tmp = Image.new("L", (int(tw)+2*pad, int(th)+2*pad), 0)
            ImageDraw.Draw(tmp).text((pad, pad + a), text, fill=255, font=font, anchor="ls")
            tmp = tmp.rotate(angle, expand=1, resample=Image.NEAREST)
            x, y = xy_from_norm(xp, yp); bx, by = tmp.size[0]//2, tmp.size[1]//2
            img.paste(tmp, (x - bx, y - by), tmp)

    def draw_with_superscript(base_txt, super_txt, xp, yp, tail_txt="",
                              angle=0, anchor="mm",
                              base_size=font_size, super_scale=0.72, super_raise=0.55, gap_px=0,
                              base_font=None, super_font=None, tail_font=None):
        """
        Compose (base)(superscript)(tail) with proper metrics so nothing clips.
        All runs use STIX Two Math (or fallbacks).
        """
        base_font = base_font or load_font_main(int(round(base_size)))
        super_font = super_font or load_font_main(int(round(base_size*super_scale)))
        tail_font = tail_font or base_font

        pad = max(12, int(0.45 * base_size))

        # widths
        canvas = Image.new("L", (2,2), 0)
        meas = ImageDraw.Draw(canvas)
        w_base = meas.textlength(base_txt, font=base_font) if base_txt else 0
        w_sup  = meas.textlength(super_txt, font=super_font) if super_txt else 0
        w_tail = meas.textlength(tail_txt, font=tail_font) if tail_txt else 0

        a_b, d_b = base_font.getmetrics()
        a_s, d_s = super_font.getmetrics()

        # vertical sizing so the superscript never clips
        raise_px = int(round(base_size * super_raise))
        top_over = max(0, (raise_px + a_s) - a_b)         # extra above baseline
        height = a_b + d_b + top_over

        Wtmp = int(math.ceil(w_base + (gap_px if super_txt else 0) + w_sup + w_tail)) + 2*pad
        Htmp = height + 2*pad
        tmp = Image.new("L", (max(1, Wtmp), max(1, Htmp)), 0)
        tdr = ImageDraw.Draw(tmp)

        # baselines
        y_base = pad + top_over + a_b
        x = pad
        if base_txt:
            tdr.text((x, y_base), base_txt, fill=255, font=base_font, anchor="ls")
            x += w_base
        if super_txt:
            tdr.text((x + gap_px, y_base - raise_px), super_txt, fill=255, font=super_font, anchor="ls")
            x += gap_px + w_sup
        if tail_txt:
            tdr.text((x, y_base), tail_txt, fill=255, font=tail_font, anchor="ls")

        if angle != 0:
            tmp = tmp.rotate(angle, expand=1, resample=Image.NEAREST)

        X, Y = xy_from_norm(xp, yp); bx, by = tmp.size[0]//2, tmp.size[1]//2
        img.paste(tmp, (X - bx, Y - by), tmp)

    # ---------------- Labels ----------------
    # Regions
    draw_text("Region I (exterior)",    0.67,  0.00)
    draw_text("Region III (exterior)", -0.67,  0.00)
    draw_text("Region II (black hole interior, r < 2M)",  0.00,  0.67)
    draw_text("Region IV (white hole interior, r < 2M)",  0.00, -0.67)

    # Singularities
    draw_text("Spacelike singularity  r = 0", 0.00,  1.0 - SING_MARGIN_FRAC + 0.03)
    draw_text("Spacelike singularity  r = 0", 0.00, -1.0 + SING_MARGIN_FRAC - 0.03)

    # Horizons / bifurcation sphere
    draw_text("Event horizons  r = 2M", 0.00, 0.00)
    draw_text("Bifurcation 2-sphere",   0.00, 0.08)

    # Rotated horizon labels H^±
    draw_with_superscript("H", "+",  0.38,  0.38, angle=-45)
    draw_with_superscript("H", "−", -0.38, -0.38, angle=-45)
    draw_with_superscript("H", "+", -0.38,  0.38, angle=+45)
    draw_with_superscript("H", "−",  0.38, -0.38, angle=+45)

    # ?^± with tail text
    draw_with_superscript("?", "+",  0.58,  0.46, tail_txt=" (future null infinity)", angle=-45)
    draw_with_superscript("?", "−",  0.58, -0.46, tail_txt=" (past null infinity)",   angle=+45)
    draw_with_superscript("?", "+", -0.58,  0.46, angle=+45)
    draw_with_superscript("?", "−", -0.58, -0.46, angle=-45)

    # Timelike/spacelike infinities: i^±, i^0
    draw_with_superscript("i", "+",  0.00,  1.02)
    draw_with_superscript("i", "−",  0.00, -1.02)
    draw_with_superscript("i", "0",  1.02,  0.00)
    draw_with_superscript("i", "0", -1.02,  0.00)

    # Asymptotically flat directions
    draw_text("r → ∞",  0.90, 0.00)
    draw_text("r → ∞", -0.90, 0.00)

    # ER bridge — note the en dash “–”
    draw_text("t = 0  (Einstein–Rosen bridge)", 0.00, -0.06)

    # Build masks
    fill_mask = (np.array(img, dtype=np.uint8) > 0)

    # Stroke = dilation of glyphs minus fill
    rad = max(1, int(round(stroke_px)))
    size = 2*rad + 1
    dilated = img.filter(ImageFilter.MaxFilter(size))
    stroke_mask = (np.array(dilated, dtype=np.uint8) > 0) & (~fill_mask)

    return fill_mask, stroke_mask

# ----------------------------
# Kruskal relations & compactification
# ----------------------------
def C_of_r(r):   # (1 - r/2) e^{r/2}
    return (1.0 - 0.5*r) * math.exp(0.5*r)

def compactify_UV_to_xy(U, V, k=KRUSKAL_K):
    Uc = (2.0/np.pi) * np.arctan(k * U)
    Vc = (2.0/np.pi) * np.arctan(k * V)
    xprime = 0.5*(Vc - Uc)
    yprime = 0.5*(Vc + Uc)
    return xprime, yprime

# ----------------------------
# EVEN spacing helpers for r-curves
# ----------------------------
def r_from_xprime_exterior(xp_target, k=KRUSKAL_K):
    # Exterior (I,III): ensure each r>2 curve hits y'=0 at evenly spaced x'.
    X_target = math.tan(0.5*math.pi * xp_target) / k
    S = X_target * X_target
    lo = 2.0 + 1e-8; hi = 6.0
    def g(r): return C_of_r(r) + S  # want C(r) = -S
    while g(hi) > 0.0:
        hi *= 1.7
        if hi > 1e6: raise RuntimeError("r_from_xprime_exterior bracket fail")
    for _ in range(80):
        mid = 0.5*(lo+hi)
        if g(mid) > 0: lo = mid
        else:          hi = mid
    return 0.5*(lo+hi)

def xprime_at_y_for_r(r, y_ref, k=KRUSKAL_K):
    """
    For r in (0,2), find x'(r) at the compactified horizontal level y'=y_ref (|y_ref|<1):
      arctan(kU) + arctan(kV) = y_ref * π, with UV = C(r) > 0.
    """
    C = C_of_r(r)
    A = y_ref * math.pi
    TA = math.tan(A)
    S = (TA / k) * (1.0 - (k*k) * C)
    disc = S*S - 4.0*C
    if disc <= 0:
        return None
    rt = math.sqrt(disc)
    U1 = 0.5*(S + rt)
    U2 = 0.5*(S - rt)
    def pick(U):
        V = S - U
        # y_ref>0 (II): U,V < 0 ; y_ref<0 (IV): U,V > 0
        sgn = -1.0 if y_ref > 0 else +1.0
        if (U*sgn > 0) and (V*sgn > 0):
            Uc = (2.0/np.pi) * math.atan(k * U)
            Vc = (2.0/np.pi) * math.atan(k * V)
            return 0.5*(Vc - Uc)
        return None
    x1 = pick(U1); x2 = pick(U2)
    if x1 is None and x2 is None:
        return None
    return x1 if x1 is not None else x2

def interior_r_list_even_x(n_int, y_ref_mag=Y_REF_INT, k=KRUSKAL_K, scan_pts=1200, edge_trim=0.06):
    """
    Build r<2 list by scanning x'(r; +y_ref_mag) over r∈(0,2), then invert by interpolation
    to place n_int curves at evenly spaced x' across the achievable band.
    """
    y_ref = float(np.clip(y_ref_mag, 0.05, 0.90))
    rmin, rmax = 1e-7, 2.0 - 1e-7
    rs = np.linspace(rmin, rmax, scan_pts)
    xs = np.array([xprime_at_y_for_r(float(r), +y_ref, k) for r in rs], dtype=np.float64)
    ok = ~np.isnan(xs)
    rs, xs = rs[ok], xs[ok]
    if xs.size < 16:
        # Fallback: evenly spaced by C in (0,2)
        Cs = np.linspace(0.95, 0.05, n_int)
        return [invert_C_to_r(float(C)) for C in Cs]
    order = np.argsort(xs)
    xs = xs[order]; rs = rs[order]
    x_lo, x_hi = xs[0], xs[-1]
    span = x_hi - x_lo
    if span < 0.2:
        Cs = np.linspace(0.95, 0.05, n_int)
        return [invert_C_to_r(float(C)) for C in Cs]
    x_lo2 = x_lo + edge_trim * span
    x_hi2 = x_hi - edge_trim * span
    x_targets = np.linspace(x_lo2, x_hi2, n_int)
    r_targets = np.interp(x_targets, xs, rs)
    for i in range(1, len(r_targets)):
        if not (r_targets[i] > r_targets[i-1] + 1e-8):
            r_targets[i] = r_targets[i-1] + 1e-6
    return list(map(float, r_targets))

def invert_C_to_r(C):
    # Monotone on (0,2): bisection
    lo, hi = 1e-10, 2.0 - 1e-10
    def f(r): return C_of_r(r) - C
    for _ in range(80):
        mid = 0.5*(lo+hi)
        if f(mid) > 0: lo = mid
        else:          hi = mid
    return 0.5*(lo+hi)

def choose_even_r_lists(n_ext, n_int, k=KRUSKAL_K, edge_xp=0.06, y_ref_int=Y_REF_INT):
    # Exterior r>2 via even x' at y'=0
    x_targets_ext = np.linspace(edge_xp, 1.0 - edge_xp, n_ext, dtype=np.float64)
    r_ext = [r_from_xprime_exterior(float(xp), k=k) for xp in x_targets_ext]
    # Interior r<2 via robust monotonic inversion at mid-wedge
    r_int = interior_r_list_even_x(n_int, y_ref_mag=y_ref_int, k=k)
    return r_ext, r_int

# ----------------------------
# Constant-t list (even “angle”)
# ----------------------------
def build_even_t_list(n_t, theta_max_deg=THETA_MAX_DEG):
    if n_t <= 0: return []
    th_max = math.radians(theta_max_deg)
    thetas = np.linspace(-th_max, th_max, n_t)
    out = []
    for th in thetas:
        m = math.tan(th)
        if abs(m) >= 1.0:
            m = math.copysign(0.999999, m)
        t = 4.0 * np.arctanh(m)
        out.append(float(t))
    return out

# ----------------------------
# Grid rasterization (sub-pixel horizon margin)
# ----------------------------
def make_coord_grid_mask(masks, line_px, n_samples,
                         r_ext_list, r_int_list, t_list, k_scale=KRUSKAL_K,
                         eps_hor_uc=1e-6, eps_inf_uc=0.003):
    """
    eps_hor_uc: tiny exclusion around U=0 in compactified Uc units (sub-pixel).
    eps_inf_uc: exclusion near |Uc|=1 to avoid degeneracy near corners.
    """
    xp_grid, yp_grid, half, cx, cy = masks["xp"], masks["yp"], masks["half"], masks["cx"], masks["cy"]
    H, W = xp_grid.shape

    img = PILImage.new("L", (W, H), 0)
    draw = ImageDraw.Draw(img)

    def to_px(xp, yp):
        x = np.rint(cx + xp * half).astype(np.int32)
        y = np.rint(cy - yp * half).astype(np.int32)
        return x, y

    def draw_parametric_from_Uc(Uc_lo, Uc_hi, V_from_U, keep_region=None):
        # Chebyshev nodes → denser near ends
        t = np.linspace(0.0, 1.0, n_samples)
        Uc = 0.5*((Uc_hi+Uc_lo) + (Uc_hi-Uc_lo)*np.cos((1.0 - t)*math.pi))
        U = (1.0/k_scale) * np.tan(0.5*np.pi * Uc)
        V = V_from_U(U)
        if keep_region:
            uv = U * V
            ok = (uv < 0.0) if keep_region == "ext" else (uv > 0.0)
            U = U[ok]; V = V[ok]
        finite = np.isfinite(U) & np.isfinite(V)
        U, V = U[finite], V[finite]
        if U.size < 2: return
        xp, yp = compactify_UV_to_xy(U, V, k=k_scale)
        keep = (np.abs(xp) + np.abs(yp) <= 1.02)
        if np.count_nonzero(keep) < 2: return
        xpx, ypx = to_px(xp[keep].astype(np.float32), yp[keep].astype(np.float32))
        pts = list(map(tuple, np.stack([xpx, ypx], axis=1)))
        if len(pts) >= 2:
            draw.line(pts, fill=255, width=line_px, joint="curve")

    # Constant-r hyperbolae: UV = C(r)
    for r in r_ext_list:  # I and III (C<0)
        C = C_of_r(r)
        draw_parametric_from_Uc(-1.0+eps_inf_uc, -eps_hor_uc, lambda U: C / U)  # Region I
        draw_parametric_from_Uc( eps_hor_uc,  1.0-eps_inf_uc, lambda U: C / U)  # Region III

    for r in r_int_list:  # II and IV (C>0)
        C = C_of_r(r)
        draw_parametric_from_Uc(-1.0+eps_inf_uc, -eps_hor_uc, lambda U: C / U)  # Region II
        draw_parametric_from_Uc( eps_hor_uc,  1.0-eps_inf_uc, lambda U: C / U)  # Region IV

    # Constant-t (region-aware)
    for t in t_list:
        eta = math.exp(0.5*t)
        draw_parametric_from_Uc(-1.0+eps_inf_uc, -eps_hor_uc, lambda U: -eta*U, keep_region="ext")  # I
        draw_parametric_from_Uc( eps_hor_uc,  1.0-eps_inf_uc, lambda U: -eta*U, keep_region="ext")  # III
        draw_parametric_from_Uc(-1.0+eps_inf_uc, -eps_hor_uc, lambda U:  eta*U, keep_region="int")  # II
        draw_parametric_from_Uc( eps_hor_uc,  1.0-eps_inf_uc, lambda U:  eta*U, keep_region="int")  # IV

    return (np.array(img, dtype=np.uint8) > 0) & masks["diamond"]

# ----------------------------
# Optional: sparse curved non-radial null rays
# ----------------------------
def r_star_of_r(r):  return r + 2.0 * np.log(abs(r/2.0 - 1.0))
def f_null_radial_factor(b, r): return 1.0 - (1.0 - 2.0/r) * (b*b) / (r*r)
def dt_dr(b, r):
    f = f_null_radial_factor(b, r)
    return 1.0 / ((1.0 - 2.0/r) * np.sqrt(np.maximum(f, 1e-16)))

def find_turning_radius(b, r_hi=2000.0):
    lo, hi = 3.0, r_hi
    def f(r): return f_null_radial_factor(b, r)
    if f(hi) <= 0: raise RuntimeError("Failed to bracket turning radius")
    if f(lo) > 0: lo = 3.0001
    for _ in range(80):
        mid = 0.5*(lo + hi)
        if f(mid) <= 0: lo = mid
        else:           hi = mid
    return 0.5*(lo+hi)

def compactify_uv_to_xy(u, v, k=COMPACT_K_NULL):
    Uc = (2.0/np.pi) * np.arctan(k * u)
    Vc = (2.0/np.pi) * np.arctan(k * v)
    xprime = 0.5*(Vc - Uc); yprime = 0.5*(Vc + Uc)
    return xprime, yprime

def make_curved_null_mask(masks, width_px, b_list, rmax, samples, k_scale):
    if not b_list:
        return np.zeros_like(masks["diamond"], dtype=bool)
    xp_grid, yp_grid, half, cx, cy = masks["xp"], masks["yp"], masks["half"], masks["cx"], masks["cy"]
    H, W = xp_grid.shape
    img = PILImage.new("L", (W, H), 0)
    draw = ImageDraw.Draw(img)

    def to_px(xp, yp):
        x = np.rint(cx + xp * half).astype(np.int32)
        y = np.rint(cy - yp * half).astype(np.int32)
        return x, y

    for b in b_list:
        r_turn = find_turning_radius(b, r_hi=max(rmax, 100.0))
        y = np.linspace(0.0, 1.0, samples, dtype=np.float64)
        r_half = r_turn + (rmax - r_turn) * (y**2)
        integrand = dt_dr(b, r_half)
        t_half = np.cumsum(0.5*(integrand[1:] + integrand[:-1]) * (r_half[1:] - r_half[:-1]))
        t_half = np.concatenate(([0.0], t_half))
        rstar_half = r_star_of_r(r_half)
        rstar_in  = rstar_half[::-1];  t_in  = -t_half[::-1]
        rstar_out = rstar_half;        t_out =  t_half
        rstar = np.concatenate([rstar_in, rstar_out[1:]])
        t     = np.concatenate([t_in,     t_out[1:]])
        u = t - rstar; v = t + rstar
        xp, yp = compactify_uv_to_xy(u, v, k=k_scale)
        keep = (np.abs(xp) + np.abs(yp) <= 1.02)
        xp = xp[keep]; yp = yp[keep]
        if xp.size < 2: continue
        xpx, ypx = to_px(xp.astype(np.float32), yp.astype(np.float32))
        pts = list(map(tuple, np.stack([xpx, ypx], axis=1)))
        draw.line(pts, fill=255, width=width_px, joint="curve")

    return (np.array(img, dtype=np.uint8) > 0) & masks["diamond"]

# ----------------------------
# Region gradients (lightest at region centers)
# ----------------------------
def region_focal_points():
    return {
        "I":   ( +0.50,  0.00),
        "III": ( -0.50,  0.00),
        "II":  (  0.00, +0.50),
        "IV":  (  0.00, -0.50),
    }

def lab_L_for_rgb(cols):
    rgb01 = cols.astype(np.float32) / 255.0
    L, _, _ = rgb01_to_lab(rgb01)
    return L

def assign_region_gradient(flat, idxs, cols, xp, yp, focus_xy):
    if len(idxs) == 0: return
    h, w = yp.shape
    yy, xx = np.divmod(idxs, w)
    xps = xp[yy, xx]; yps = yp[yy, xx]
    fx, fy = focus_xy
    d = np.hypot(xps - fx, yps - fy)
    wts = 1.0 - np.clip(d, 0.0, 1.2) / 1.2
    wts = wts**1.2
    order_pix = np.argsort(-wts)          # closest first
    Ls = lab_L_for_rgb(cols)
    order_col = np.argsort(-Ls)           # brightest first
    flat[idxs[order_pix]] = cols[order_col]

# ----------------------------
# Background (outside diamond): low-chroma ripples
# ----------------------------
def assign_background_pattern(flat, idxs, cols, xp, yp):
    if len(idxs) == 0: return
    h, w = yp.shape
    yy, xx = np.divmod(idxs, w)
    xps = xp[yy, xx]; yps = yp[yy, xx]
    r1 = np.abs(xps) + np.abs(yps)                # diamond radius
    ang = np.arctan2(yps, xps)
    key = r1 + 0.020*np.sin(18.0*r1 + 4.0*ang)    # subtle ripple
    order_pix = np.argsort(key)
    rgb01 = cols.astype(np.float32) / 255.0
    L, a, b = rgb01_to_lab(rgb01)
    C = np.sqrt(a*a + b*b)
    order_col = np.lexsort((np.abs(L-55.0), C))   # prefer lower chroma, mid L
    flat[idxs[order_pix]] = cols[order_col]

# ----------------------------
# Fill helper
# ----------------------------
def fill_indices(img_flat, idxs, colors):
    assert colors.shape[0] == idxs.shape[0]
    img_flat[idxs] = colors

# ----------------------------
# Main
# ----------------------------
def main():
    t0 = time.time()

    masks = make_masks(W, H, MARGIN_FRAC, LINE_PX, SING_MARGIN_FRAC)
    xp, yp, half = masks["xp"], masks["yp"], masks["half"]

    # Base masks (exclusivity comes later)
    label_fill_mask, label_stroke_mask = place_labels_masks(W, H, masks, LABEL_FONT_SIZE, LABEL_STROKE_PIX)
    line_mask_base = (masks["horizons"] | masks["singular"])  # within diamond already
    frame_band = (np.abs(np.abs(xp) + np.abs(yp) - 1.0) <= (FRAME_PX / half))  # diamond outline band

    # Grid / t-lines
    if DRAW_GRID:
        r_ext_list, r_int_list = choose_even_r_lists(N_R_EXT, N_R_INT, k=KRUSKAL_K, y_ref_int=Y_REF_INT)
        t_list = build_even_t_list(N_T, THETA_MAX_DEG)
        grid_mask_raw = make_coord_grid_mask(
            masks, GRID_LINE_PX, GRID_SAMPLES,
            r_ext_list, r_int_list, t_list,
            k_scale=KRUSKAL_K,
            eps_hor_uc=1e-6,    # sub-pixel near horizons
            eps_inf_uc=0.003    # tight near corners
        )
    else:
        grid_mask_raw = np.zeros_like(masks["diamond"], dtype=bool)

    # Optional non-radial null geodesics
    if DRAW_CURVED_NULL:
        geo_mask_raw = make_curved_null_mask(
            masks, width_px=GEO_CURVE_PX,
            b_list=GEODESIC_IMPACT_PARAMS,
            rmax=GEODESIC_RMAX,
            samples=GEODESIC_SAMPLES,
            k_scale=COMPACT_K_NULL
        )
    else:
        geo_mask_raw = np.zeros_like(masks["diamond"], dtype=bool)

    # Region base fills and outside
    reg_I_base, reg_II_base = masks["I"], masks["II"]
    reg_III_base, reg_IV_base = masks["III"], masks["IV"]
    bg_base = masks["outside"]

    # -------- Exclusive partition (TOP → DOWN) --------
    occupied = np.zeros((H, W), dtype=bool)

    def take(mask):
        m = mask & (~occupied)
        occupied[:] = occupied | m
        return m

    # Visual precedence: labels → frame → lines → grid → geo → regions → background
    m_label_fill = take(label_fill_mask)
    m_label_stk  = take(label_stroke_mask)
    m_frame      = take(frame_band)
    m_line       = take(line_mask_base)
    m_grid       = take(grid_mask_raw)
    m_geo        = take(geo_mask_raw)
    m_I          = take(reg_I_base)
    m_II         = take(reg_II_base)
    m_III        = take(reg_III_base)
    m_IV         = take(reg_IV_base)
    m_bg         = take(bg_base)

    assert occupied.all(), "Partition left some pixels unassigned"

    # Pixel counts
    n_label_fill = int(np.count_nonzero(m_label_fill))
    n_label_stk  = int(np.count_nonzero(m_label_stk))
    n_frame      = int(np.count_nonzero(m_frame))
    n_line       = int(np.count_nonzero(m_line))
    n_grid       = int(np.count_nonzero(m_grid))
    n_geo        = int(np.count_nonzero(m_geo))
    n_I          = int(np.count_nonzero(m_I))
    n_II         = int(np.count_nonzero(m_II))
    n_III        = int(np.count_nonzero(m_III))
    n_IV         = int(np.count_nonzero(m_IV))
    n_bg         = int(np.count_nonzero(m_bg))

    total = (n_label_fill + n_label_stk + n_frame + n_line + n_grid + n_geo
             + n_I + n_II + n_III + n_IV + n_bg)
    assert total == W*H, f"Mask coverage error: {total} vs {W*H}"

    print("Pixel budget (exclusive):")
    print(f"  Label fill   : {n_label_fill:,} (lightest)")
    print(f"  Label stroke : {n_label_stk:,} (darkest)")
    print(f"  Frame        : {n_frame:,} (diamond outline)")
    print(f"  Lines        : {n_line:,} (horizons + singularities)")
    print(f"  Grid         : {n_grid:,} (even r- and t-lines)")
    print(f"  Geo          : {n_geo:,} (optional non-radial null)")
    print(f"  Regions I-IV : {n_I:,}, {n_II:,}, {n_III:,}, {n_IV:,}")
    print(f"  Background   : {n_bg:,}")
    print(f"  TOTAL        : {total:,}")

    # Palette
    pal = BinPalette()

    # LIGHTEST → label fill
    cols_label_fill = pal.allocate_light(n_label_fill)
    # DARKEST → label outline, frame, horizons/singularities, grid, optional null
    cols_label_stk  = pal.allocate_dark(n_label_stk)
    cols_frame      = pal.allocate_dark(n_frame)
    cols_lines      = pal.allocate_dark(n_line)
    cols_grid       = pal.allocate_dark(n_grid)
    cols_geo        = pal.allocate_dark(n_geo)
    # Low-chroma → background (patterned)
    cols_bg = pal.allocate_low_chroma(n_bg)

    # Region colors (sorted by hue preferences)
    bins_I   = pal.bins_by_hue_preference(REGION_HUE_PREFS["I"],   MIN_CHROMA_FOR_REGIONS)
    bins_III = pal.bins_by_hue_preference(REGION_HUE_PREFS["III"], MIN_CHROMA_FOR_REGIONS)
    bins_II  = pal.bins_by_hue_preference(REGION_HUE_PREFS["II"],  MIN_CHROMA_FOR_REGIONS)
    bins_IV  = pal.bins_by_hue_preference(REGION_HUE_PREFS["IV"],  MIN_CHROMA_FOR_REGIONS)
    cols_I   = pal.allocate_from_bins(bins_I,   n_I)
    cols_III = pal.allocate_from_bins(bins_III, n_III)
    cols_II  = pal.allocate_from_bins(bins_II,  n_II)
    cols_IV  = pal.allocate_from_bins(bins_IV,  n_IV)

    # Assemble
    img = np.empty((H, W, 3), dtype=np.uint8)
    flat = img.reshape(-1, 3)

    def idx(mask): return np.flatnonzero(mask.ravel())

    idx_label_fill = idx(m_label_fill)
    idx_label_stk  = idx(m_label_stk)
    idx_frame      = idx(m_frame)
    idx_line       = idx(m_line)
    idx_grid       = idx(m_grid)
    idx_geo        = idx(m_geo)
    idx_I          = idx(m_I)
    idx_II         = idx(m_II)
    idx_III        = idx(m_III)
    idx_IV         = idx(m_IV)
    idx_bg         = idx(m_bg)

    # Background (outside diamond): patterned, subdued
    assign_background_pattern(flat, idx_bg, cols_bg, xp, yp)

    # Region gradients (lightest at region centers)
    foc = region_focal_points()
    assign_region_gradient(flat, idx_I,   cols_I,   xp, yp, foc["I"])
    assign_region_gradient(flat, idx_III, cols_III, xp, yp, foc["III"])
    assign_region_gradient(flat, idx_II,  cols_II,  xp, yp, foc["II"])
    assign_region_gradient(flat, idx_IV,  cols_IV,  xp, yp, foc["IV"])

    # Overlays (exclusive masks, order doesn't change totals)
    fill_indices(flat, idx_grid,        cols_grid)
    fill_indices(flat, idx_geo,         cols_geo)
    fill_indices(flat, idx_line,        cols_lines)
    fill_indices(flat, idx_frame,       cols_frame)
    fill_indices(flat, idx_label_stk,   cols_label_stk)
    fill_indices(flat, idx_label_fill,  cols_label_fill)

    if VERIFY_UNIQUENESS:
        packed = (img[...,0].astype(np.uint32) << 16) | (img[...,1].astype(np.uint32) << 8) | img[...,2].astype(np.uint32)
        uniq = np.unique(packed)
        assert uniq.size == 256*256*256, f"Uniqueness failed: {uniq.size} unique colors"
        print("Verified: all 16,777,216 colors used exactly once.")

    ts = datetime.now().strftime("%Y%m%d_%H%M%S")
    out_path = os.path.join(OUTPUT_DIR, f"penrose_allrgb_{W}x{H}_{ts}.png")
    PILImage.fromarray(img, mode="RGB").save(out_path, compress_level=0)
    print(f"Saved: {out_path}  (elapsed {time.time() - t0:.1f}s)")

if __name__ == "__main__":
    main()