# Colab-ready: AllRGB (4096x4096) encoding of the EHT/ESO polarized M87* image (eso2105a)
# Bias: dark pixels -> neutral near-black colors; bright pixels -> orange-ish colors
# Still uses every 24-bit RGB value exactly once.
import os
import gc
import requests
import numpy as np
from PIL import Image
# -----------------------------
# Settings
# -----------------------------
OUT_SIZE = 4096
N = OUT_SIZE * OUT_SIZE # 16,777,216
ESO_TIF = "https://cdn.eso.org/images/original/eso2105a.tif"
ESO_JPG = "https://cdn.eso.org/images/large/eso2105a.jpg"
DOWNLOAD_PATH = "eso2105a_source"
OUT_PNG = "m87_polarized_allrgb_orange_black_4096.png"
# 0..65280 (because we use integer Rec.709 weights summing to 256)
# Increase to push "neutral blacks" deeper into midtones; decrease to make orange bias kick in earlier.
DARK_LUMA_CUTOFF = 12000
# -----------------------------
# Helpers
# -----------------------------
def download(url: str, path: str, chunk=1<<20):
print("Downloading:", url)
r = requests.get(url, stream=True, timeout=60)
r.raise_for_status()
with open(path, "wb") as f:
for b in r.iter_content(chunk_size=chunk):
if b:
f.write(b)
return path
def load_source_image():
tif_path = DOWNLOAD_PATH + ".tif"
jpg_path = DOWNLOAD_PATH + ".jpg"
if not os.path.exists(tif_path) and not os.path.exists(jpg_path):
try:
download(ESO_TIF, tif_path)
path = tif_path
except Exception as e:
print("TIFF download failed, falling back to JPEG.\nReason:", repr(e))
download(ESO_JPG, jpg_path)
path = jpg_path
else:
path = tif_path if os.path.exists(tif_path) else jpg_path
img = Image.open(path).convert("RGB")
print("Loaded:", path, "size:", img.size)
return img
def center_crop_to_square(img: Image.Image) -> Image.Image:
w, h = img.size
s = min(w, h)
left = (w - s) // 2
top = (h - s) // 2
return img.crop((left, top, left + s, top + s))
def resize_to_4096(img: Image.Image) -> Image.Image:
return img.resize((OUT_SIZE, OUT_SIZE), resample=Image.LANCZOS)
def luma16_from_rgb_u8(rgb_u8: np.ndarray) -> np.ndarray:
"""
Luma in 0..65280 using integer weights summing to 256:
luma16 = 54*R + 183*G + 19*B
"""
r = rgb_u8[..., 0].astype(np.uint32)
g = rgb_u8[..., 1].astype(np.uint32)
b = rgb_u8[..., 2].astype(np.uint32)
return (54*r + 183*g + 19*b).astype(np.uint16)
def xor_u32_range(n: int) -> int:
m = n - 1
def xor_0_to(m):
r = m & 3
if r == 0: return m
if r == 1: return 1
if r == 2: return m + 1
return 0
return xor_0_to(m)
def orange_score_int(r, g, b):
"""
Higher = more orange-ish (warm):
prefer B low, R high, G mid, and R~G.
Works on uint32 arrays; returns int16-ish range.
"""
r = r.astype(np.int32)
g = g.astype(np.int32)
b = b.astype(np.int32)
score = (2*r + g) - (4*b) - np.abs(r - g)
# extra bias for "R>=G>=B"
score += 256 * (r >= g).astype(np.int32)
score += 256 * (g >= b).astype(np.int32)
return score.astype(np.int32)
def neutral_score_int(r, g, b):
"""
Higher = more neutral/gray (channels close).
"""
r = r.astype(np.int32)
g = g.astype(np.int32)
b = b.astype(np.int32)
d = np.abs(r - g) + np.abs(g - b) # 0..510
return (510 - d).astype(np.int32)
# -----------------------------
# 1) Load and prep the scientific image
# -----------------------------
src = load_source_image()
src_sq = center_crop_to_square(src)
src_4096 = resize_to_4096(src_sq)
target_rgb = np.asarray(src_4096, dtype=np.uint8) # (4096,4096,3)
target_luma16 = luma16_from_rgb_u8(target_rgb) # uint16 0..65280
# Compute per-pixel style key:
# - dark pixels: prefer neutralness
# - bright pixels: prefer orange-ness
Rpx = target_rgb[..., 0].astype(np.uint32)
Gpx = target_rgb[..., 1].astype(np.uint32)
Bpx = target_rgb[..., 2].astype(np.uint32)
px_orange = orange_score_int(Rpx, Gpx, Bpx)
px_neutral = neutral_score_int(Rpx, Gpx, Bpx)
# Style score switches based on darkness of the *pixel*:
px_style = np.where(target_luma16 < DARK_LUMA_CUTOFF, px_neutral, px_orange).astype(np.int32)
# Free PIL objects
del src, src_sq, src_4096
gc.collect()
# -----------------------------
# 2) Pixel order: (luma ascending, style descending within same luma)
# -----------------------------
print("Sorting pixels by (luma, style)...")
# lexsort: last key is primary -> primary=luma, secondary=-style
pix_order = np.lexsort((-px_style.reshape(-1), target_luma16.reshape(-1)))
del Rpx, Gpx, Bpx, px_orange, px_neutral, px_style
gc.collect()
# -----------------------------
# 3) Palette: allRGB values 0..2^24-1, sorted by (luma, style)
# Dark luma -> prefer neutral blacks; bright luma -> prefer orange-ish
# -----------------------------
print("Building allRGB values 0..2^24-1 ...")
vals = np.arange(N, dtype=np.uint32)
r = (vals >> 16).astype(np.uint32)
g = ((vals >> 8) & 255).astype(np.uint32)
b = (vals & 255).astype(np.uint32)
pal_luma16 = (54*r + 183*g + 19*b).astype(np.uint16)
col_orange = orange_score_int(r, g, b)
col_neutral = neutral_score_int(r, g, b)
col_style = np.where(pal_luma16 < DARK_LUMA_CUTOFF, col_neutral, col_orange).astype(np.int32)
print("Sorting palette by (luma, style)...")
col_order = np.lexsort((-col_style, pal_luma16))
colors_sorted = vals[col_order] # uint32 length N
del r, g, b, pal_luma16, col_orange, col_neutral, col_style, col_order, vals
gc.collect()
# -----------------------------
# 4) Assign: darkest pixels -> darkest colors, with the style bias within each luma
# -----------------------------
print("Assigning colors to pixels...")
out_flat = np.empty(N, dtype=np.uint32)
out_flat[pix_order] = colors_sorted
del pix_order, colors_sorted, target_rgb, target_luma16
gc.collect()
# -----------------------------
# 5) Unpack to uint8 RGB image and save
# -----------------------------
print("Unpacking to RGB and saving PNG...")
R = ((out_flat >> 16) & 255).astype(np.uint8)
G = ((out_flat >> 8) & 255).astype(np.uint8)
B = (out_flat & 255).astype(np.uint8)
out_rgb = np.stack([R, G, B], axis=1).reshape(OUT_SIZE, OUT_SIZE, 3)
# Fast permutation sanity checks (not a full uniqueness scan)
expected_sum = (np.uint64(N - 1) * np.uint64(N)) // 2
actual_sum = out_flat.astype(np.uint64).sum()
assert actual_sum == expected_sum, "Sum check failed (colors likely not a permutation)."
expected_xor = np.uint32(xor_u32_range(N))
actual_xor = np.bitwise_xor.reduce(out_flat.astype(np.uint32))
assert actual_xor == expected_xor, "XOR check failed (colors likely not a permutation)."
img_out = Image.fromarray(out_rgb, mode="RGB")
img_out.save(OUT_PNG, compress_level=3)
print("Wrote:", OUT_PNG)
# Optional preview (Colab)
try:
import matplotlib.pyplot as plt
plt.figure(figsize=(6, 6))
plt.imshow(out_rgb)
plt.axis("off")
plt.show()
except Exception:
pass
# Optional: download from Colab
# from google.colab import files
# files.download(OUT_PNG)
