#!/usr/bin/env python3
"""The Jupiter aperture ladder, photographed rather than simulated.

Eight rungs, eight instruments, eight pipelines: five amateur telescopes from
150 to 408 mm, Hubble at 2.4 m above the atmosphere, Gemini North at 8.1 m
using the very lucky-imaging-and-stacking technique the book calls "how to
manipulate a photo", and Juno at the planet. Every disc is found automatically
(largest bright blob), cropped and resampled so Jupiter's equatorial diameter
is the same number of pixels in every tile; nothing is sharpened here. The
Juno rung is a close pass and is not to the common scale.

The source images are third parties' and are not kept in this repository;
they are used with each photographer's permission and credited on the figure.
Point --src at a folder holding them under the names below.

    python3 build_jupiter_ladder.py --src /path/to/_ladder  -> docs/planetary-imaging/img/fig-jupiter-ladder.png
"""
import argparse
import os

import numpy as np
from PIL import Image, ImageDraw, ImageFilter, ImageFont
from scipy import ndimage

HERE = os.path.dirname(os.path.abspath(__file__))
OUT = os.path.join(HERE, "..", "docs", "planetary-imaging", "img", "fig-jupiter-ladder.png")
GO = os.path.join(HERE, "..", "docs", "planetary-imaging", "img", "fig-go-stack.png")

# rung: (file, aperture label, instrument line, date/technique line, credit line)
RUNGS = [
    ("07_ASK_kowollik_150mm_j2013-10-25.jpg", "150 mm", "Maksutov-Cassegrain, Ludwigsburg", "25 Oct 2013 · stacked, sharpened", "Silvia Kowollik"),
    ("05_ASK_peralta_152mm_j2015-01-24_GRS_visible.png", "152 mm", "Meade LS 6″ ACF, Algeciras", "24 Jan 2015 · AutoStakkert!", "Javier Peralta"),
    ("06_ASK_shanos_250mm_j2026-06-24.jpg", "250 mm", "Meade LX200, Sarasota · 685 nm IR", "24 Jun 2026 · one frame, not stacked", "Gregory T. Shanos"),
    ("GO", "356 mm", "Celestron C14, Cebu", "22 Jul 2022 · best 5% of 3,215 frames", "Christopher Go, raw capture; stacked for this page"),
    ("08_ASK_barry_406mm_j2016-07-17_GRS_visible.png", "408 mm", "Newtonian, Broken Hill · 685 nm IR", "17 Jul 2016 · stacked, sharpened", "Trevor Barry"),
    ("01_hubble_2.4m_heic1410a.jpg", "2.4 m · Hubble", "WFC3, in orbit", "21 Apr 2014 · no seeing to beat", "NASA / ESA / A. Simon (GSFC), CC BY 4.0"),
    ("02_gemini_8.1m_noirlab2011a.jpg", "8.1 m · Gemini North", "NIRI, Maunakea · 4.7 µm", "2019 · lucky imaging, best 10% stacked", "Int'l Gemini Obs. / NOIRLab / NSF / AURA, M. H. Wong et al."),
    ("04_juno_processed_thomopoulos_id11625.jpg", "Juno · at the planet", "JunoCam, perijove 21", "21 Jul 2019 · km per pixel — not to scale", "NASA / JPL-Caltech / SwRI / MSSS · processing Thomopoulos, CC BY"),
]
DISC_PX = 300      # equatorial diameter in every tile
TILE = 400
PAD = 18


# the Gemini frame is thermal infrared: the disc is bright only in its belts, so it needs a
# lower threshold and a wider blur to be seen as one blob
THRESH = {"02_gemini_8.1m_noirlab2011a.jpg": (0.05, 6)}


def find_disc(im, thr_frac=0.18, blur=2):
    """Bounding box of the planet: the largest bright blob that is round -- a bounding box
    near square and a fill fraction near pi/4 -- so caption panels and text are passed over."""
    g = np.asarray(im.convert("L").filter(ImageFilter.GaussianBlur(blur)), dtype=float)
    thr = thr_frac * g.max()
    lab, n = ndimage.label(g > thr)
    best = None
    for k in range(1, n + 1):
        ys, xs = np.where(lab == k)
        w, h = xs.max() - xs.min() + 1, ys.max() - ys.min() + 1
        if w < 40 or h < 40:
            continue
        aspect, fill = w / h, len(xs) / (w * h)
        if 0.8 < aspect < 1.25 and 0.6 < fill < 0.92:
            if best is None or len(xs) > best[0]:
                best = (len(xs), xs.min(), ys.min(), xs.max(), ys.max())
    if best is None:
        raise SystemExit("no round bright blob found")
    return best[1:]


def mask_outside(tile, radius_px):
    """Black out everything beyond about 1.1 disc radii, with a soft edge: captions, compass roses, frame borders."""
    w, h = tile.size
    yy, xx = np.mgrid[0:h, 0:w]
    r = np.hypot(xx - w / 2 + 0.5, yy - h / 2 + 0.5)
    m = np.clip((1.12 * radius_px - r) / (0.05 * radius_px), 0, 1)
    a = np.asarray(tile, dtype=float) * m[:, :, None]
    return Image.fromarray(a.astype(np.uint8))


def go_tile():
    """The best-5% stack from fig-go-stack.png, middle panel."""
    im = Image.open(GO).convert("RGB")
    w, h = im.size
    return im.crop((int(w * 0.352), int(h * 0.158), int(w * 0.652), int(h * 0.86)))


def tile_for(path, label):
    if path == "GO":
        im = go_tile()
    else:
        im = Image.open(path).convert("RGB")
    if label.startswith("Juno"):
        w, h = im.size
        c = im.crop((int(w * 0.42), int(h * 0.0), int(w * 0.96), int(h * 0.54)))
        return c.resize((TILE, TILE), Image.LANCZOS)
    x0, y0, x1, y1 = find_disc(im, *THRESH.get(os.path.basename(path), (0.18, 2)))
    cx, cy = (x0 + x1) / 2, (y0 + y1) / 2
    d = x1 - x0                      # equatorial width in source pixels
    s = DISC_PX / d
    half = TILE / 2 / s
    box = (int(cx - half), int(cy - half), int(cx + half), int(cy + half))
    canvas = Image.new("RGB", im.size, "black")
    canvas.paste(im, (0, 0))
    # pad so the crop never leaves the image
    big = Image.new("RGB", (im.size[0] + 2 * int(half), im.size[1] + 2 * int(half)), "black")
    big.paste(canvas, (int(half), int(half)))
    box = tuple(int(v + half) for v in box)
    return mask_outside(big.crop(box).resize((TILE, TILE), Image.LANCZOS), DISC_PX / 2)


def font(sz, bold=False):
    for f in (["/usr/share/fonts/truetype/dejavu/DejaVuSans-Bold.ttf"] if bold else ["/usr/share/fonts/truetype/dejavu/DejaVuSans.ttf"]):
        if os.path.exists(f):
            return ImageFont.truetype(f, sz)
    return ImageFont.load_default()


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--src", required=True)
    ap.add_argument("--out", default=OUT)
    a = ap.parse_args()
    cols, rows = 4, 2
    cap_h = 92
    W = cols * (TILE + PAD) + PAD
    H = 74 + rows * (TILE + cap_h + PAD) + PAD
    sheet = Image.new("RGB", (W, H), "white")
    dr = ImageDraw.Draw(sheet)
    dr.text((PAD, 16), "The same planet up the ladder — eight instruments, eight pipelines", fill="#181715", font=font(26, True))
    dr.text((PAD, 48), f"Each image as its author released it; the only step here is resampling every disc to {DISC_PX} px across the equator. No sharpening or processing has been added to any of them.",
            fill="#6E675B", font=font(15))
    for i, (fn, lab, inst, when, cred) in enumerate(RUNGS):
        r, c = divmod(i, cols)
        x = PAD + c * (TILE + PAD)
        y = 74 + r * (TILE + cap_h + PAD)
        t = tile_for(fn if fn == "GO" else os.path.join(a.src, fn), lab)
        sheet.paste(t, (x, y))
        dr.rectangle((x, y, x + TILE - 1, y + TILE - 1), outline="#bdb6a7")
        dr.text((x, y + TILE + 8), lab, fill="#2E4057", font=font(19, True))
        dr.text((x, y + TILE + 34), inst, fill="#3A352D", font=font(13))
        dr.text((x, y + TILE + 52), when, fill="#3A352D", font=font(13))
        dr.text((x, y + TILE + 70), cred, fill="#6E675B", font=font(11))
    os.makedirs(os.path.dirname(a.out), exist_ok=True)
    sheet.save(a.out, optimize=True)
    print("wrote", os.path.normpath(a.out), sheet.size)


if __name__ == "__main__":
    main()