#!/usr/bin/env python3
# Builds wwwroot/vendor/globe/earth-blue-marble.jpg from NASA's Blue Marble Next Generation (topography and bathymetry,
# December 2004, public domain), downloaded once:
#   https://eoimages.gsfc.nasa.gov/images/imagerecords/73000/73909/world.topo.bathy.200412.3x5400x2700.jpg
# usage: tools/globe/prepare-day-texture.py <that jpg> <output jpg>      (needs Pillow; not part of the build)
#
# NASA's mosaic does not wrap cleanly: towards the poles a strip at its east edge is lighter than the map on both sides
# of the 180° meridian, which the globe draws as a line from each pole. On every row where the map agrees on both sides
# of the seam but the strip around it does not, the strip is replaced by a blend between the two sides. A real feature
# on the meridian (Fiji, the Ross Ice Shelf) makes the two sides differ, so its rows are left alone. Within half a degree
# of each pole a row is a circle a few kilometres across, so it gets one colour, its median, which removes the strip
# there too. Then the map is scaled to 4096×2048 (the globe's texture size) with Lanczos.
import sys
from PIL import Image

BAND = 16      # pixels on each side of the seam that may be repaired
ANCHOR = 4     # pixels averaged just outside the band, on each side
AGREE = 24     # the two sides agree when their mean channels differ by less than this
ODD = 30       # a band pixel is an artefact when it differs from the blend by more than this
POLAR = 0.5    # degrees from a pole within which a row is one colour


def mean(pixels):
    return tuple(sum(p[i] for p in pixels) / len(pixels) for i in range(3))


def distance(a, b):
    return sum(abs(x - y) for x, y in zip(a, b)) / 3


source, target = sys.argv[1], sys.argv[2]
image = Image.open(source).convert("RGB")
width, height = image.size
pixels = image.load()
repaired = 0
for y in range(height):
    west = mean([pixels[width - BAND - 1 - i, y] for i in range(ANCHOR)])
    east = mean([pixels[BAND + i, y] for i in range(ANCHOR)])
    if distance(west, east) >= AGREE:
        continue
    columns = [width - BAND + i for i in range(BAND)] + list(range(BAND))
    blend = [tuple(w + (e - w) * (i + 1) / (2 * BAND + 1) for w, e in zip(west, east)) for i in range(2 * BAND)]
    if max(distance(pixels[x, y], b) for x, b in zip(columns, blend)) <= ODD:
        continue
    for x, b in zip(columns, blend):
        pixels[x, y] = tuple(round(c) for c in b)
    repaired += 1

polar = 0
for y in range(height):
    latitude = 90 - (y + 0.5) * 180 / height
    if abs(latitude) < 90 - POLAR:
        continue
    row = [pixels[x, y] for x in range(width)]
    median = tuple(sorted(p[i] for p in row)[len(row) // 2] for i in range(3))
    for x in range(width):
        pixels[x, y] = median
    polar += 1

image.resize((4096, 2048), Image.LANCZOS).save(target, "JPEG", quality=90, optimize=True, progressive=True)
print(f"{repaired} of {height} rows repaired at the seam, {polar} polar rows made one colour; wrote {target}")