sky-cam/moon_composite.py

#!/usr/bin/env python3
"""moon_composite.py — paint a high-res lunar texture into a sky-cam frame.

Goal: from a frame east captured of a small white blob (~38 px), produce a
1920x1080 image that looks like east took it through a 65x telephoto.

What's real, from east:
  - Sky color, atmospheric halo, any clouds drifting past
  - Time, parallactic angle (orientation of "up" on the moon)
  - Position of the moon in the frame at that instant

What's borrowed:
  - The lunar surface texture (one cached high-res reference image)

Library entry point:

  from moon_composite import composite_full_moon
  composite_full_moon(
      source_jpg='/data/east/2026-04-29/21-07-00.jpg',
      detection=detect_moon(...),
      when_utc=datetime(2026, 4, 29, 21, 7, 0, tzinfo=timezone.utc),
      ref_moon_path='/path/to/full-moon.jpg',
      out_path='/movies/east/full-moons/2026-04-full-moon.jpg',
  )

CLI:

  python3 moon_composite.py SOURCE.jpg WHEN_UTC REF_MOON.jpg OUT.jpg
"""
from __future__ import annotations

import argparse
import math
import sys
from datetime import datetime, timezone
from pathlib import Path

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

# Pillow ≥ 10 renamed resampling constants
try:
    LANCZOS = Image.Resampling.LANCZOS
    BICUBIC = Image.Resampling.BICUBIC
except AttributeError:  # Pillow < 10
    LANCZOS = Image.LANCZOS
    BICUBIC = Image.BICUBIC

# Common DejaVu Sans paths across distros; first hit wins.
_DEJAVU_CANDIDATES = [
    '/usr/share/fonts/truetype/dejavu/DejaVuSans.ttf',
    '/usr/share/fonts/dejavu/DejaVuSans.ttf',
    '/usr/share/fonts/TTF/DejaVuSans.ttf',
    '/Library/Fonts/DejaVuSans.ttf',  # macOS homebrew
]


def _load_font(size: int = 18) -> ImageFont.ImageFont:
    """Return a Unicode-capable font at `size` pt, or the PIL default."""
    for path in _DEJAVU_CANDIDATES:
        try:
            return ImageFont.truetype(path, size)
        except (OSError, IOError):
            pass
    return ImageFont.load_default()


def _draw_caption(draw: ImageDraw.ImageDraw, caption: str,
                  xy: tuple[int, int], output_w: int) -> None:
    """Draw a drop-shadow caption, sanitising characters the default font
    cannot encode when no TTF is available."""
    font = _load_font(max(14, output_w // 100))
    # If we fell back to the bitmap default font it only handles latin-1; strip
    # anything outside that range rather than crashing.
    if not hasattr(font, 'path'):
        caption = caption.encode('latin-1', errors='replace').decode('latin-1')
        caption = caption.replace('\x3f', '-')  # '?' placeholder → dash
    x, y = xy
    draw.text((x + 2, y + 2), caption, fill=(0, 0, 0), font=font)
    draw.text((x,     y),     caption, fill=(220, 220, 220), font=font)


def _square_crop_to_disk(ref: Image.Image, threshold: int = 25) -> Image.Image:
    """Tight-crop a reference moon image to the disk's bounding square.

    Most lunar reference photos have the disk on a large black field with
    significant padding.  We threshold on luminance and crop to bbox + small
    margin so the disk fills our target square evenly.
    """
    g = np.asarray(ref.convert('L'))
    mask = g >= threshold
    ys, xs = np.where(mask)
    if len(xs) == 0:
        return ref
    x0, x1 = int(xs.min()), int(xs.max())
    y0, y1 = int(ys.min()), int(ys.max())
    cx = (x0 + x1) // 2
    cy = (y0 + y1) // 2
    half = max(x1 - x0, y1 - y0) // 2 + 4  # tiny margin
    left = max(0, cx - half)
    top = max(0, cy - half)
    right = min(ref.width, cx + half)
    bottom = min(ref.height, cy + half)
    return ref.crop((left, top, right, bottom))


def _disk_mask(size: int, feather_px: int = 6) -> Image.Image:
    """Soft circular alpha mask the size of the reference moon image."""
    m = Image.new('L', (size, size), 0)
    draw = ImageDraw.Draw(m)
    # Inset slightly so the feather sits inside the disk edge
    draw.ellipse((feather_px, feather_px, size - feather_px, size - feather_px),
                 fill=255)
    if feather_px > 0:
        m = m.filter(ImageFilter.GaussianBlur(radius=feather_px))
    return m


def _apply_phase_shadow(disk: Image.Image, phase_angle_deg: float,
                        waxing: bool) -> Image.Image:
    """Darken the un-lit portion of the disk based on phase.

    phase_angle_deg: 0 = full, 90 = quarter, 180 = new.
    waxing: True = lit on right, False = lit on left.

    Implementation: the terminator is an ellipse whose semi-minor axis is
    cos(phase_angle).  Pixels on the un-lit side are multiplied by a small
    factor (not zero, so the un-lit limb stays visible like real earthshine).
    """
    if phase_angle_deg < 1.0:
        return disk  # full enough that shadow would be a single-pixel sliver
    w, h = disk.size
    # Build a mask: 1.0 in lit area, 0.04 in un-lit area, soft transition near
    # the terminator.
    cx, cy = w / 2.0, h / 2.0
    r = min(w, h) / 2.0
    yy, xx = np.mgrid[0:h, 0:w].astype(np.float32)
    # Normalise to disk coords (-1..1)
    nx = (xx - cx) / r
    ny = (yy - cy) / r
    # Distance from disk centre (we still want to clip to the disk)
    in_disk = (nx * nx + ny * ny) <= 1.0
    # Terminator equation: x_norm = cos(phase)  on the appropriate side.
    # For waxing moon, the lit portion is right of the terminator (nx > x_t).
    cos_p = math.cos(math.radians(phase_angle_deg))
    # When the moon is more than half lit (cos_p > 0), terminator is on the
    # un-lit side and the lit portion is broader.  When less than half
    # (cos_p < 0), terminator is on the lit side.
    # Distance from terminator (positive = lit side)
    if waxing:
        d = nx - (-cos_p)
    else:
        d = -(nx - cos_p)
    # Smooth step around the terminator (~2% of radius)
    soft_px = max(1.5 / r, 0.01)
    lit = np.clip(0.5 + d / (2 * soft_px), 0.04, 1.0)
    lit = np.where(in_disk, lit, 1.0)  # leave outside-disk untouched
    arr = np.asarray(disk).astype(np.float32)
    arr = arr * lit[..., None]
    return Image.fromarray(np.clip(arr, 0, 255).astype(np.uint8), disk.mode)


def composite_full_moon(
    source_jpg: str,
    detection,                                  # MoonDetection from moon_detect
    when_utc: datetime,
    ref_moon_path: str,
    out_path: str,
    output_size: tuple[int, int] = (1920, 1080),
    moon_height_pct: float = 0.70,
    caption: str | None = None,
):
    """Build the full-moon close-up composite and write it to out_path."""
    # Lazy import — avoids loading skyfield when caller doesn't need it
    import moon_phase

    src = Image.open(source_jpg).convert('RGB')
    cx, cy = detection.centroid_xy
    src_diam = detection.diameter_px

    # ── Crop east around the moon, sized so the moon fills moon_height_pct ──
    crop_h = int(round(src_diam / moon_height_pct))
    crop_w = int(round(crop_h * output_size[0] / output_size[1]))
    sw, sh = src.size
    crop_w = min(crop_w, sw)
    crop_h = min(crop_h, sh)
    left = int(round(cx - crop_w / 2))
    top = int(round(cy - crop_h / 2))
    left = max(0, min(left, sw - crop_w))
    top = max(0, min(top, sh - crop_h))
    crop = src.crop((left, top, left + crop_w, top + crop_h))
    bg = crop.resize(output_size, LANCZOS)

    # ── Where is the moon's center within the upscaled background? ──
    moon_x_in_crop = cx - left
    moon_y_in_crop = cy - top
    scale = output_size[1] / crop_h
    out_moon_cx = moon_x_in_crop * scale
    out_moon_cy = moon_y_in_crop * scale

    # ── Load reference texture, tight-crop to disk ──
    ref = Image.open(ref_moon_path).convert('RGB')
    ref = _square_crop_to_disk(ref)
    target_size = int(round(output_size[1] * moon_height_pct))
    target_size += target_size % 2  # even
    ref_resized = ref.resize((target_size, target_size), LANCZOS)

    # ── Phase shadow (skip when essentially full) ──
    illum = moon_phase.illumination(when_utc)
    pa = moon_phase.phase_angle(when_utc)
    if illum < 0.995:
        wax = moon_phase.waxing(when_utc)
        ref_resized = _apply_phase_shadow(ref_resized, pa, wax)

    # ── Parallactic-angle rotation ──
    par = moon_phase.parallactic_angle(when_utc)
    # PIL rotates counter-clockwise for positive angles; we want celestial
    # north to end up "up" in the camera image.  Negate so the rotation
    # direction matches image-space y-down convention.
    ref_rot = ref_resized.rotate(-par, resample=BICUBIC, expand=False)

    # ── Composite with feathered circular mask ──
    feather = max(4, target_size // 200)
    mask = _disk_mask(target_size, feather_px=feather)

    paste_x = int(round(out_moon_cx - target_size / 2))
    paste_y = int(round(out_moon_cy - target_size / 2))
    # Clamp so the disk stays fully on canvas (recenter if needed)
    paste_x = max(0, min(paste_x, output_size[0] - target_size))
    paste_y = max(0, min(paste_y, output_size[1] - target_size))
    bg.paste(ref_rot, (paste_x, paste_y), mask)

    # ── Caption ──
    if caption:
        draw = ImageDraw.Draw(bg)
        _draw_caption(draw, caption, (22, output_size[1] - 48), output_size[0])

    Path(out_path).parent.mkdir(parents=True, exist_ok=True)
    bg.save(out_path, quality=92)
    return out_path


def _make_atmosphere_layer(
    east_frame_path: str,
    output_size: tuple[int, int],
    blur_radius: int = 0,
) -> Image.Image:
    """Scale the full east frame to output_size and blur to atmospheric haze.

    The blur removes wide-angle camera detail (RTSP artefacts, OSD text,
    pixel noise) while preserving real sky colour and large-scale cloud
    structure.  The resulting layer is applied OVER the NASA moon disk so
    it reads as "clouds between the observer and the moon" — which is
    physically correct.

    blur_radius=0 → auto: output_width // 10, which smooths pixel-level
    detail but keeps cloud-scale gradients visible.
    """
    src = Image.open(east_frame_path).convert('RGB')
    layer = src.resize(output_size, LANCZOS)
    r = blur_radius if blur_radius > 0 else output_size[0] // 10
    return layer.filter(ImageFilter.GaussianBlur(radius=r))


def render_phase_closeup(
    nasa_render_path: str,
    out_path: str,
    output_size: tuple[int, int] = (1920, 1080),
    moon_height_pct: float = 0.92,
    caption: str | None = None,
    background: tuple[int, int, int] = (0, 0, 0),
    east_frame_path: str | None = None,
    atmosphere_opacity: float = 0.0,
    atmosphere_blur: int = 0,
):
    """Full-screen close-up rendering using a NASA SVS Dial-a-Moon image.

    Rendering pipeline:
      1. Black background (outer space).
      2. NASA moon disk — correct phase, libration, crater shadows — centred
         and scaled to moon_height_pct of frame height.
      3. Atmospheric layer (optional): east's camera frame, scaled to output
         size and blurred, composited OVER the moon at atmosphere_opacity.
         This is physically correct — clouds are between the observer and the
         moon so they occlude the disk, not sit behind it.

    atmosphere_opacity controls what the viewer sees through east's sky:
      0.00 — perfectly clear: pure NASA render, no overlay
      0.10 — slight haze: moon is sharp but slightly softened
      0.35 — noticeable cloud cover: moon partially obscured
      0.65 — heavy overcast: moon a faint glow through thick cloud
    """
    # ── Background + NASA moon ────────────────────────────────────────────
    bg = Image.new('RGB', output_size, background)
    moon = Image.open(nasa_render_path).convert('RGB')
    moon = _square_crop_to_disk(moon)

    target = int(round(output_size[1] * moon_height_pct))
    target += target % 2
    moon_resized = moon.resize((target, target), LANCZOS)

    feather = max(3, target // 240)
    mask = _disk_mask(target, feather_px=feather)

    px = (output_size[0] - target) // 2
    py = (output_size[1] - target) // 2
    bg.paste(moon_resized, (px, py), mask)

    # ── Atmospheric layer from east frame — applied OVER the moon ─────────
    if east_frame_path and atmosphere_opacity > 0.0:
        atm = _make_atmosphere_layer(east_frame_path, output_size, atmosphere_blur)
        bg = Image.blend(bg, atm, alpha=atmosphere_opacity)

    if caption:
        draw = ImageDraw.Draw(bg)
        _draw_caption(draw, caption, (22, output_size[1] - 48), output_size[0])

    Path(out_path).parent.mkdir(parents=True, exist_ok=True)
    bg.save(out_path, quality=92)
    return out_path


def _cli():
    p = argparse.ArgumentParser()
    p.add_argument('source')
    p.add_argument('when_utc', help='ISO 8601 UTC, e.g. 2026-04-29T21:07:00Z')
    p.add_argument('ref_moon')
    p.add_argument('out')
    p.add_argument('--caption', default=None)
    p.add_argument('--moon-pct', type=float, default=0.70)
    args = p.parse_args()

    from moon_detect import detect_moon
    det = detect_moon(args.source)
    if det is None:
        print('ERROR: no moon detected in source', file=sys.stderr)
        return 2
    when = datetime.fromisoformat(args.when_utc.replace('Z', '+00:00'))
    composite_full_moon(
        args.source, det, when, args.ref_moon, args.out,
        moon_height_pct=args.moon_pct, caption=args.caption,
    )
    print(f'wrote {args.out}')
    return 0


if __name__ == '__main__':
    sys.exit(_cli())