#!/usr/bin/env python3
"""
M87* all-RGB image (4096×4096), scientifically-informed crescent.
What this does
--------------
1) Builds a physically-motivated intensity map I(x,y):
- Thin Gaussian ring ("photon ring" scale) with a strong SOUTHERN brightness
asymmetry (relativistic beaming) and a deep central brightness depression.
- Parameters guided by EHT measurements (ring diameter ~42 μas, thin ring).
Refs: EHT 2019 reports; PRIMO 2023 sharper ring; 2024 multi-epoch ring
with southerly brightness. See citations in the chat message.
2) Uses ALL 24-bit RGB colors exactly once:
- Total pixels = 4096×4096 = 16,777,216.
- Selects exactly as many “ring pixels” as needed by a geometric ring mask.
- Assigns the most chromatic colors to the ring (ordered by hue along angle).
- Assigns the remaining colors to all other pixels, dark-to-bright to match I.
Design choices
--------------
* Ring geometry:
- r0_frac ~ 0.29 (radius as fraction of image size), thin width.
- Southern bright sector via a Doppler-like factor A(θ) = (1+k cos(θ-θ0))^p,
with θ0 = south (−π/2) in image coordinates (y down).
* Exact-use guarantee:
- Colors are split into two non-overlapping pools (ring vs non-ring).
- Each pool is strictly permutation-assigned to its pixel set.
- A final uniqueness assertion checks we used every RGB exactly once.
Runtime & memory
----------------
This is heavy by design (tens to hundreds of MB of RAM and big sorts).
For lower RAM, you can chunk some steps, but this single-pass version
is easiest to verify for exact color accounting.
Outputs
-------
- PNG named: M87_allrgb_.png
- Console prints: counts & verification.
"""
import numpy as np
from PIL import Image
import time
import math
# ------------------------- parameters you might tweak -------------------------
SIZE = 4096 # image width=height; must be 4096 to match 16,777,216 colors
SAVE_TIMESTAMP = time.strftime("%Y%m%d_%H%M%S")
OUTFILE = f"M87_allrgb_{SAVE_TIMESTAMP}.png"
# Ring geometry (fraction of SIZE)
r0_frac = 0.29 # ring radius
ring_band = 96.0 # half-thickness of the ring band in pixels (visual ring thickness)
sigma_px = 28.0 # Gaussian radial sigma controlling the ring profile thickness
# Beaming (brightness asymmetry, brightest toward south)
theta0 = -math.pi/2 # south (in image coords: +y down)
k_beam = 0.65 # 0..1 strength of cosine beaming
p_beam = 2.2 # raises contrast of beaming
# Background/contrast shaping
bg_envelope_strength = 0.25 # faint background glow
shadow_sharpness = 8.0 # higher -> steeper dark center
# -----------------------------------------------------------------------------
# --- 1) Build coordinate grid and physically-motivated intensity I(x,y) ---
H = W = SIZE
yy, xx = np.indices((H, W), dtype=np.float32)
cx = (W - 1) * 0.5
cy = (H - 1) * 0.5
# Pixel-space radius & angle (θ measured from +x to -y for "image up")
dx = xx - cx
dy = yy - cy
r = np.hypot(dx, dy)
theta = np.arctan2(-dy, dx) # so that -pi/2 points downward (south) in the image
r0 = r0_frac * SIZE
# Radial "thin ring": Gaussian around r0
I_rad = np.exp(-0.5 * ((r - r0) / sigma_px) ** 2)
# Doppler-like beaming: brightest near theta0 (south)
beam = 1.0 + k_beam * np.cos(theta - theta0)
beam = np.clip(beam, 0.0, None) ** p_beam
# Deep central depression (shadow) via a sharp logistic around r ~ r0 - a few sigmas
# This suppresses the interior brightness by >10:1 relative to ring.
inner_edge = r0 - 2.2 * sigma_px
shadow = 1.0 / (1.0 + np.exp(-(r - inner_edge) / (sigma_px / shadow_sharpness)))
# Faint background envelope to avoid a totally flat outer region
bg = bg_envelope_strength * np.exp(- (r / (0.42 * SIZE)) ** 2)
# Combine and normalize
I = (I_rad * beam * shadow) + bg
I -= I.min()
I /= (I.max() + 1e-12)
# --- ring mask: geometric band around r0 (uses band, not intensity threshold) ---
ring_mask = (np.abs(r - r0) <= ring_band) & (I > 0.02)
ring_idx = np.flatnonzero(ring_mask)
non_idx = np.flatnonzero(~ring_mask)
N_total = SIZE * SIZE
N_ring = ring_idx.size
N_non = non_idx.size
assert N_ring + N_non == N_total
# Order ring pixels around azimuth, starting at theta0 (south) and going CCW
phi = (theta - theta0) % (2.0 * np.pi)
phi_ring = phi.flat[ring_idx]
order_ring_pixels = ring_idx[np.argsort(phi_ring, kind='mergesort')]
# Order non-ring pixels by intensity (dark -> bright)
I_non = I.flat[non_idx]
order_non_pixels = non_idx[np.argsort(I_non, kind='mergesort')]
# --- 2) Generate ALL 24-bit colors (no repeats), then split into two pools ---
# Colors as uint8 channels via bit tricks; stable and memory-light(ish)
N = 1 << 24 # 16,777,216
ids = np.arange(N, dtype=np.uint32)
R = ((ids >> 16) & 0xFF).astype(np.uint8)
G = ((ids >> 8) & 0xFF).astype(np.uint8)
B = ( ids & 0xFF).astype(np.uint8)
# "Chromaticity" metric to pick the most colorful colors for the ring:
# Using max(R,G,B) - min(R,G,B) (channel spread). This equals (S*V)*255 in sRGB scale.
maxc = np.maximum(np.maximum(R, G), B).astype(np.int16)
minc = np.minimum(np.minimum(R, G), B).astype(np.int16)
spread = (maxc - minc).astype(np.int16) # 0..255
# Select EXACTLY N_ring "most chromatic" colors for the ring
# (Tie-breaking handled by argpartition+argsort for determinism.)
if N_ring > 0:
part_idx = np.argpartition(spread, -N_ring)[-N_ring:]
# Now order those ring colors by hue for a smooth hue cycle around the ring
# Compute hue ONLY for the selected subset to save memory/time.
r_f = R[part_idx].astype(np.float32) / 255.0
g_f = G[part_idx].astype(np.float32) / 255.0
b_f = B[part_idx].astype(np.float32) / 255.0
mx = np.maximum(np.maximum(r_f, g_f), b_f)
mn = np.minimum(np.minimum(r_f, g_f), b_f)
df = mx - mn
# Hue in [0,1)
h = np.zeros_like(mx, dtype=np.float32)
mask = df > 0
# Avoid division by zero; compute hue sectors
r_is_max = (mx == r_f) & mask
g_is_max = (mx == g_f) & mask
b_is_max = (mx == b_f) & mask
# Note: vectorized hue (0..6) then /6
h[r_is_max] = ( (g_f[r_is_max] - b_f[r_is_max]) / df[r_is_max] ) % 6.0
h[g_is_max] = ( (b_f[g_is_max] - r_f[g_is_max]) / df[g_is_max] ) + 2.0
h[b_is_max] = ( (r_f[b_is_max] - g_f[b_is_max]) / df[b_is_max] ) + 4.0
h *= (1.0 / 6.0)
order_ring_colors_local = np.argsort(h, kind='mergesort')
ring_color_ids = part_idx[order_ring_colors_local]
else:
ring_color_ids = np.array([], dtype=np.uint32)
# The remaining colors go to non-ring pixels, ordered by luminance (dark -> bright).
# Integer luma approximation: Y ≈ 0.2126R + 0.7152G + 0.0722B, scaled by 10,000 for integers.
if N_non > 0:
used_mask = np.zeros(N, dtype=bool)
used_mask[ring_color_ids] = True
non_color_ids = np.flatnonzero(~used_mask)
# Luma as integer to avoid huge float arrays
Yint = (2126 * R[non_color_ids].astype(np.uint32) +
7152 * G[non_color_ids].astype(np.uint32) +
722 * B[non_color_ids].astype(np.uint32)) # 0..2,550,000
order_non_local = np.argsort(Yint, kind='mergesort')
nonring_color_ids = non_color_ids[order_non_local]
else:
nonring_color_ids = np.array([], dtype=np.uint32)
# --- 3) Assemble the final image buffer (bijective mapping) ---
img = np.empty((N_total, 3), dtype=np.uint8)
# Assign ring: pixels ordered by azimuth, colors ordered by hue
if N_ring > 0:
Rc = (ring_color_ids >> 16) & 0xFF
Gc = (ring_color_ids >> 8) & 0xFF
Bc = (ring_color_ids ) & 0xFF
img[order_ring_pixels, 0] = Rc.astype(np.uint8)
img[order_ring_pixels, 1] = Gc.astype(np.uint8)
img[order_ring_pixels, 2] = Bc.astype(np.uint8)
# Assign non-ring: pixels ordered by intensity (dark->bright), colors ordered by luma
if N_non > 0:
Rn = (nonring_color_ids >> 16) & 0xFF
Gn = (nonring_color_ids >> 8) & 0xFF
Bn = (nonring_color_ids ) & 0xFF
img[order_non_pixels, 0] = Rn.astype(np.uint8)
img[order_non_pixels, 1] = Gn.astype(np.uint8)
img[order_non_pixels, 2] = Bn.astype(np.uint8)
# Reshape and save
img = img.reshape((H, W, 3))
Image.fromarray(img, mode="RGB").save(OUTFILE, compress_level=0)
# --- 4) Verify exact usage of all 16,777,216 colors ---
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 == (1 << 24), f"Color count mismatch: {uniq.size}"
print("Success!")
print(f"Ring pixels: {N_ring:,} | Non-ring pixels: {N_non:,}")
print(f"Saved: {OUTFILE}")
Why this maps to the real M87*
- The ring with a southern bright sector and deep central depression mirrors what EHT reported (ring diameter ≈ 42 μas, strong azimuthal asymmetry consistent with relativistic beaming). The 2023 PRIMO reconstruction tightened the ring width; we emulate a thin band. The 2024 multi-epoch analysis again found a southerly brightness distribution and a ring size consistent with 42 μas.
