#!/usr/bin/env python3
"""The photographs on the long-path page, annotated with what was measured on them.

  1. A frame of the book's own Chicago time-lapse (Nowicki, 30 April 2015, 16 s in):
     the water horizon, the level the four roofs share, the Willis tiers and masts,
     and the cut the fit puts on the city.
  2. The wide frame of the Pontchartrain clip the book links (25 s in): the ten
     resolved towers, the straight line through their tops, the pile caps.
  3. The p. 111 still, cropped to the skyline: the overlay's 1,014-ft line as printed,
     against the hidden heights a beach camera actually has.

The source frames are third parties' and are reproduced for critical review; they
are not kept in this repository. Point --src at a folder holding the two videos and
the extracted p. 111 image under the names below.

    python3 annotate_long_path_frames.py --src /path/to/folder
"""
import argparse
import math
import os
import subprocess
import sys
import tempfile

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import chicago_refraction as cr      # noqa: E402  the frame fit
import long_path_k_bands as kb       # noqa: E402  the hiding formula

from PIL import Image, ImageDraw, ImageEnhance, ImageFont

HERE = os.path.dirname(os.path.abspath(__file__))
IMG = os.path.join(HERE, "..", "docs", "long-path-cases", "img")
NOWICKI = "chicago_michigan_timelapse_raw.mkv"
PONTCH = "pontchartrain_raw.mkv"
STILL = "p111-still.png"          # the 800x457 image embedded at p. 111

INK, GOLD, RED, PALE = "#ffffff", "#f2c14e", "#ff6b57", "#dfe6ee"


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


def frame_at(video, second, td):
    """The same one-frame-per-second sampling the measurement scripts use (fps=1), so the
    frame numbers here are the frame numbers there."""
    d = os.path.join(td, os.path.basename(video) + "_frames"); os.makedirs(d, exist_ok=True)
    if not os.listdir(d):
        subprocess.run(["ffmpeg", "-loglevel", "error", "-y", "-i", video, "-vf", "fps=1", os.path.join(d, "f%03d.png")], check=True)
    return Image.open(os.path.join(d, f"f{second + 1:03d}.png")).convert("RGB")


def label(d, xy, text, fill, sz=15, bold=False, anchor="la", box=True):
    f = font(sz, bold)
    if box:
        x0, y0, x1, y1 = d.textbbox(xy, text, font=f, anchor=anchor)
        d.rectangle((x0 - 4, y0 - 2, x1 + 4, y1 + 2), fill=(0, 0, 0, 140))
    d.text(xy, text, fill=fill, font=f, anchor=anchor)


def nowicki(src, td):
    im = frame_at(os.path.join(src, NOWICKI), 16, td)
    wide = im.resize((960, 540), Image.LANCZOS)
    wide.save(os.path.join(IMG, "fig-nowicki-wide.jpg"), quality=86)
    print("wrote fig-nowicki-wide.jpg", wide.size)
    im = ImageEnhance.Contrast(im).enhance(1.8)
    crop = (300, 560, 1560, 900)                      # the skyline band
    c = im.crop(crop)
    c = c.resize((c.width, c.height), Image.LANCZOS)
    ov = Image.new("RGBA", c.size, (0, 0, 0, 0)); d = ImageDraw.Draw(ov)
    hz = 825 - crop[1]; roof = hz - 88; wroof = hz - 137; wtip = hz - 182
    sec, cut, mpp, err = next(f for f in cr.fit_frames() if f[0] == 16)
    k_here = (kb.k_globe(cr.H_OBS, cr.D_WILLIS, cut + err), kb.k_globe(cr.H_OBS, cr.D_WILLIS, cut - err))
    # where the water line would sit for other coefficients: the same buildings, more or less of them hidden
    for k, name, col in ((0.0, "airless globe, k = 0", "#9aa7b8"), (0.13, "k = 0.13, survey standard", "#9aa7b8"), (0.20, "k = 0.20, top of normal marine", "#9aa7b8")):
        C = kb.hidden_globe(cr.H_OBS, cr.D_WILLIS, k)
        yk = hz - (C - cut) / mpp
        for x in range(0, c.width, 14):
            d.line((x, yk, x + 7, yk), fill=col, width=1)
        label(d, (665, yk - 2), f"water line if {name}: {C:.0f} m hidden", col, 12, anchor="ld")
    # the band the water line swept over the usable frames of the clip
    fits = [f for f in cr.fit_frames()]
    lo, hi = min(f[1] for f in fits), max(f[1] for f in fits)
    y_lo, y_hi = hz - (hi - cut) / mpp, hz - (lo - cut) / mpp
    d.rectangle((0, y_lo, c.width, y_hi), fill=(242, 193, 78, 60))
    label(d, (8, y_hi + 4), f"over the clip (10–22 s) the water line ranged across this band: {lo:.0f}–{hi:.0f} m hidden", GOLD, 12, anchor="la")
    d.line((0, hz, c.width, hz), fill=PALE, width=2)
    label(d, (8, hz + 6), f"water horizon in this frame: {cut:.0f} ± {err:.0f} m of the city hidden, k ≈ {k_here[0]:.2f}–{k_here[1]:.2f}", PALE, 15)
    d.line((0, roof, c.width, roof), fill=GOLD, width=1)
    label(d, (c.width - 8, roof - 4), "Hancock, Aon and Trump roofs (344–357 m), Willis 90th-floor tier (360 m)", GOLD, 13, anchor="rd")
    wx = 485 - crop[0]
    d.line((wx + 45, wroof, wx + 200, wroof), fill=GOLD, width=1); label(d, (wx + 206, wroof), "Willis roof, 442 m", GOLD, 14, anchor="lm")
    d.line((wx + 30, wtip, wx + 200, wtip), fill=GOLD, width=1); label(d, (wx + 206, wtip), "Willis masts, 527 m", GOLD, 14, anchor="lm")
    for x, name in ((485, "Willis"), (725, "Aon"), (920, "Trump"), (1390, "Hancock")):
        label(d, (x - crop[0], hz + 56), name, INK, 15, bold=True, anchor="ma")
    label(d, (8, 8), "Nowicki, Grand Mere State Park, 30 April 2015, 16 s into the time-lapse the book links at p. 111 (shapedebate.com/32)", INK, 14)
    label(d, (8, 30), "camera 59 m above the lake, Willis Tower 90.8 km; scale 2.27 m per pixel across the frame", INK, 13)
    out = Image.alpha_composite(c.convert("RGBA"), ov).convert("RGB")
    out.save(os.path.join(IMG, "fig-nowicki-frame.jpg"), quality=88)
    print("wrote fig-nowicki-frame.jpg", out.size)


def horizon_fit(im):
    """The sky-to-water transition across the open water right of the towers, column band by
    column band: the row of steepest darkening, and the rows at which the darkening is a
    quarter and three-quarters done (the width of the band). Least squares through the
    steepest rows gives the horizon as a tilted line -- the camera is not level."""
    import numpy as np
    a = np.asarray(im.convert("L"), dtype=float)
    xs, ys, widths = [], [], []
    for x0 in range(700, 1240, 60):
        band = a[300:400, x0:x0 + 60].mean(1)
        g = np.diff(band)
        r = int(np.argmin(g)) + 1
        hi, lo = band[:5].mean(), band[-5:].mean()
        q1 = int(np.argmax(band < hi - 0.25 * (hi - lo))); q3 = int(np.argmax(band < hi - 0.75 * (hi - lo)))
        xs.append(x0 + 30); ys.append(300 + r); widths.append(q3 - q1)
    n = len(xs); sx = sum(xs); sy = sum(ys); sxx = sum(x * x for x in xs); sxy = sum(x * y for x, y in zip(xs, ys))
    m = (n * sxy - sx * sy) / (n * sxx - sx * sx); b = (sy - m * sx) / n
    return m, b, sum(widths) / n


def pontchartrain(src, td):
    im = frame_at(os.path.join(src, PONTCH), 25, td)
    m, b, band = horizon_fit(im)
    crop = (0, 0, 900, 480)
    c = im.crop(crop)
    ov = Image.new("RGBA", c.size, (0, 0, 0, 0)); d = ImageDraw.Draw(ov)
    # the soft band the sky-to-water transition occupies, then the fitted horizon through it
    poly = [(0, b - band / 2), (c.width, m * c.width + b - band / 2), (c.width, m * c.width + b + band / 2), (0, b + band / 2)]
    d.polygon(poly, fill=(223, 230, 238, 60))
    d.line((0, b, c.width, b + m * c.width), fill=PALE, width=2)
    tops = [(233, 101), (305, 138), (356, 164), (394, 184), (421, 199), (446, 211), (464, 221), (480, 229), (493, 237), (505, 242)]
    n = len(tops); sx = sum(x for x, _ in tops); sy = sum(y for _, y in tops); sxx = sum(x * x for x, _ in tops); sxy = sum(x * y for x, y in tops)
    mt = (n * sxy - sx * sy) / (n * sxx - sx * sx); bt = (sy - mt * sx) / n
    d.line((200, mt * 200 + bt, 640, mt * 640 + bt), fill=GOLD, width=1)
    for x, y in tops:
        d.ellipse((x - 4, y - 4, x + 4, y + 4), outline=GOLD, width=2)
    label(d, (560, 200), "the ten resolved tops: straight, as both", GOLD, 14)
    label(d, (560, 220), "models predict over the first 3 km", GOLD, 14)
    label(d, (560, 250), "the deciding towers, 10–25 km out: unresolved", RED, 14)
    label(d, (8, 400), f"water horizon (fitted): tilted {abs(math.degrees(math.atan(m))):.1f}°, the camera is not level; a soft band about {band:.0f} px deep", PALE, 14)
    label(d, (8, 422), "the far towers' bases run into that band before their tops do; nothing below it can be read", PALE, 13)
    label(d, (8, 8), "The clip the book links at p. 112 (shapedebate.com/34): Mathew Brown via dcforce, undated, 25 s in", INK, 14)
    label(d, (8, 30), "I-10 / I-55 interchange, west end of Lake Pontchartrain; the two nearest towers are a taller anchor design", INK, 13)
    out = Image.alpha_composite(c.convert("RGBA"), ov).convert("RGB")
    out.save(os.path.join(IMG, "fig-pontchartrain-frame.jpg"), quality=88)
    print("wrote fig-pontchartrain-frame.jpg", out.size, f"horizon tilt {math.degrees(math.atan(m)):.2f} deg, band {band:.0f} px")


def still(src):
    im = Image.open(os.path.join(src, STILL)).convert("RGB")
    c = im.crop((300, 250, 800, 372)); s = 2
    c = ImageEnhance.Contrast(c).enhance(1.5).resize((c.width * s, c.height * s), Image.LANCZOS)
    ov = Image.new("RGBA", c.size, (0, 0, 0, 0)); d = ImageDraw.Draw(ov)
    # the overlay's own line sits at (orig) y = 325; the breakwater cut at y = 348; the 1354 label line at 303
    for y, txt, col in ((303, "the overlay's “1,354 ft” line: the Skydeck, not the 1,451-ft roof", GOLD), (325, "the overlay's “1,014 ft” line: 8 in/mi² × 39², a surface drop, not a hidden height", RED)):
        yy = (y - 250) * s
        d.line((0, yy, c.width, yy), fill=col, width=1)
        label(d, (c.width - 8, yy - 3), txt, col, 13, anchor="rd")
    yy = (348 - 250) * s
    d.line((0, yy, c.width, yy), fill=PALE, width=2)
    label(d, (8, yy + 5), "foreground breakwater: cuts the skyline before the water horizon does", PALE, 13)
    label(d, (8, 8), "The p. 111 still, skyline enlarged 2×: a frame from a video the book does not name, with its overlay", INK, 13)
    label(d, (8, 28), "from a 2 m beach a globe hides 846 ft airless, 684 ft in standard air; the frame shows ~600–700 ft of the Willis above the breakwater", INK, 12)
    out = Image.alpha_composite(c.convert("RGBA"), ov).convert("RGB")
    out.save(os.path.join(IMG, "fig-p111-still.jpg"), quality=88)
    print("wrote fig-p111-still.jpg", out.size)


def main():
    ap = argparse.ArgumentParser()
    ap.add_argument("--src", required=True)
    a = ap.parse_args()
    os.makedirs(IMG, exist_ok=True)
    with tempfile.TemporaryDirectory() as td:
        nowicki(a.src, td)
        pontchartrain(a.src, td)
    still(a.src)


if __name__ == "__main__":
    main()