#!/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()
