M87 Black Hole v2

Unlike the first M87 Black Hole image on this site, this version uses actual eso data.

# 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)

Made with Chatgpt, Python, and Google Colab.