sky_stack_v3.py

#!/usr/bin/env -S uv run
# /// script
# requires-python = ">=3.11"
# dependencies = ["numpy", "scipy", "pillow"]
# ///
"""
sky_stack_v3.py -- build one still "decay-style" star-trail composite from a
static-camera night-sky video recording.

THE ONE IDEA THIS SCRIPT IS BUILT AROUND
-----------------------------------------
Pick ONE frame -- --base-time -- as the final image. Star trails accumulate
over the footage IMMEDIATELY BEFORE that frame, ending just short of it (the
last sample the trail accumulator actually sees sits one --sample-dt before
--base-time, not at it -- see build_star_trails()'s docstring for exactly
why):

    base   = the frame at --base-time, untouched
    trails = decay-accumulated star response over
             [base_time - trail_duration, base_time - sample_dt]
    final  = max(base, trails)

Why this is the whole trick: in a real long exposure, the bright HEAD of
every star trail is just the star itself, sitting wherever it happened to be
when the shutter closed. The base frame's own (stationary, un-trailed) stars
end up AS the trail heads -- but not because the base frame itself is fed
into the trail accumulator; it never is, it is decoded completely
separately (see grab_frame() in main()) and the trail pass never touches it.
The coincidence instead comes from the final `max(base, trails)` step:
wherever the base frame's own full-brightness, undecayed star is brighter
than whatever the trail accumulated at that pixel, the base star simply wins
outright and IS the visible trail head, and the freshest trail sample sits
close enough behind it -- at this script's recommended ~0.3px drift per
sample, see TARGET_PX_PER_SAMPLE and "SUB-PIXEL SAMPLING" below -- that the
two read as one continuous streak rather than leaving a visible gap. There
is no moment where a "raw" star sits apart from its own trail, and no
separate frame whose stars have to be matched up or erased against a
different plate later -- but that relies on --sample-dt being fine enough
for the footage's own drift rate; too coarse a step on slow-drift footage
would leave the undecayed head a full sample behind the base star instead of
seamless with it. Nothing needs cleaning up because nothing was ever wrong
to begin with.

WHAT THIS EXPECTS
------------------
  * A STATIC camera. The only thing allowed to move on screen is the
    starfield (a rotating sky dome/plane over time) and, optionally, small
    foreground subjects (a character, UI) that stay put.
  * A drifting starfield: individual stars are small (1-3px), point-like,
    and slide slowly and smoothly across the frame from sample to sample.
  * Enough LEAD-IN footage before --base-time for trails to develop -- see
    "THE ONE HARD LIMITATION" below.

If the camera pans, zooms, or cuts, none of this works: the background
percentile stack (below) assumes a given screen pixel is either "sky the
starfield sweeps across" or "static foreground" for the WHOLE accumulation
window. Stabilize the footage first if the camera moves.

THE ONE HARD LIMITATION: LEAD-IN FOOTAGE (read before picking --base-time)
----------------------------------------------------------------------------
Trails have to end exactly at the base frame, so how LONG a trail can get is
capped by how much footage exists before it, not by an independently-chosen
window elsewhere in the video. A trail needs roughly HALF_LIVES_FOR_FULL_TAIL
decay half-lives of accumulation to look fully developed (short of that, it
visibly tapers off before reaching full length, rather than fading smoothly
to nothing the way a longer trail does). Pick a --base-time too close to the
start of the footage and the trails WILL look noticeably short, no matter
how the other parameters are tuned; that is not a bug, it is exactly what
"the trail ends at the base frame" implies. The script prints how much
lead-in is actually available for the --base-time you chose, and how many
half-lives that amounts to, every time it runs -- read that line before
judging the output.

HOW STAR DETECTION WORKS
---------------------------------------------------------------------------
Every sample contributes a CONTINUOUS weight to the trail accumulator, never
a 0/1 accept-or-reject decision. A binary threshold causes two defects at
once: a star's brightness naturally "twinkles" from frame to frame in most
renderers, so a hard threshold flickers on and off around a twinkle,
punching gaps in an otherwise continuous trail; and lowering the threshold
far enough to stop losing faint stars pulls in small bright details on any
foreground subject too (hair highlights, fabric speculars), painting them
with trail noise. There is no single threshold that fixes both, because the
failure is the accept/reject step itself. So there isn't one:

  * WHITE TOP-HAT (frame minus its own grayscale morphological opening)
    keeps only small bright structures and suppresses large ones. A star
    fits inside the opening's structuring element and gets subtracted away
    to ~0; the Milky Way, clouds, a lit character, a moon disc etc. are all
    bigger than the kernel, so the opening reproduces them almost exactly
    and the top-hat difference is ~0 there too. This is the ONLY filter
    needed to reject broad structures, and it uses nothing about WHERE they
    sit on screen -- it works unchanged whether the "big bright thing" in
    frame is a galaxy band or a character's face. --tophat-k controls the
    kernel size (see its --help for how to measure it against your own
    footage's star size).
  * A soft noise floor (--thresh-lo/--thresh-hi) ramps the top-hat response
    from 0 to 1 via a smoothstep, rather than a hard cutoff, so a
    twinkling-but-real star's dim moments still contribute SOMETHING instead
    of vanishing from the trail outright, while sensor noise and compression
    artifacts sitting near-but-not-at 0 don't slowly accumulate into a grey
    haze over the whole sky.
  * CONTINUOUS ACCUMULATION: `acc = max(acc * decay, layer)` every sample,
    never a threshold-then-count. Thresholding first and only accumulating
    the survivors would make trails visibly DASHED whenever a star's
    twinkle dipped it below the cutoff for a sample or two; letting every
    sample contribute its raw (possibly small, never zero unless truly
    absent) weight and taking the max against a decaying history keeps the
    trail continuous through a twinkle's dim phase instead of punching a
    gap in it. This peak-hold shape is also why the trail head lands
    exactly on the base frame's own star (see the top of this docstring)
    and why no purely temporal filter can further reduce twinkle-driven
    modulation within one pixel's own brief transit -- see "SPATIAL
    SMOOTHING ALONG THE TRAIL DIRECTION" below for the filter that actually
    does.
  * STATIC SUPPRESSION (--bg-percentile and friends): a per-pixel LOW
    percentile of luma, taken over the whole accumulation window, is used as
    a soft multiplier on every sample's star response. The reasoning: the
    sky a starfield sweeps across is dark for MOST of the window (a star
    only occupies any given pixel briefly as it drifts through), so its low
    percentile stays close to its minimum -- near 0, i.e. "definitely sky".
    A static bright foreground detail (a character, a UI element, terrain)
    is lit for MOST of the window instead, so its low percentile stays high
    -- near its usual brightness, i.e. "definitely not sky" -- no matter how
    "star-like" (small, bright) it looks in any single frame. This is why a
    percentile and not a plain per-pixel MINIMUM: a minimum only needs ONE
    dark sample across the whole window to call a pixel "sky", and a static
    subject's idle animation (hair sway, a twitch) can supply that one dark
    sample without the pixel ever being real background even once. A low
    percentile asks "is this pixel dark MOST of the time" instead of "was it
    EVER dark", which the swaying-hair case fails and the swept-sky case
    passes. --min-static-blob-area then rejects any suppressed region too
    SMALL to plausibly be a deliberately-composed foreground subject: on a
    densely star-filled sky, a handful of unrelated stars can each pass
    close enough to the SAME open-sky pixel, at different sampled instants
    over a long window, to look statistically like one object sitting there
    continuously and tip that pixel's percentile into "suppress" -- a real
    static subject is suppressed as one large, spatially coherent blob,
    while this kind of coincidence produces small, scattered ones, so a
    single area threshold separates them without knowing anything about
    where either sits on screen (see drop_small_static_blobs() for the
    exact mechanism).
  * TRANSITION-COUNT GATE (--bg-transitions-lo/-hi): the percentile test
    above answers exactly one question -- "is this pixel dark MOST of the
    time" -- and there is a real static-camera scene where that question
    does not have a clean answer: a silhouette edge that is animated (a
    character's hair swaying in an idle animation). Such a pixel sits
    exactly on the boundary between hair and sky, so as the strand sways
    past it, the pixel alternates between "hair-dark" and "sky-dark" often
    enough that it is dark ROUGHLY HALF the time either way -- which a star
    sweeping through the same pixel also produces, briefly. On the
    "how-often-dark" axis alone these two cases are genuinely
    indistinguishable; --edge-protect-r's dilation only helps where the
    edge happens to sit next to an already-protected anchor pixel, and does
    nothing for a thin, isolated lock (a ponytail tip) with open sky on
    every side.
    The fix is a SECOND, INDEPENDENT axis: not "how often is this pixel
    dark", but "how many times does it CHANGE". A star drifting across a
    pixel is a single event -- the pixel goes dark-to-bright-to-dark once
    (maybe twice, at the sub-pixel margins) over the whole window, however
    long that window is. A swaying hair edge is a repeating event -- the
    strand sways back and forth many times over the same window, and the
    pixel flips with it every time. "Once" and "many times" are the same
    regardless of how dark the pixel happens to be on any given sample, so
    this axis separates the two populations that the percentile axis
    cannot. The two thresholds ramp between "release" and "force-suppress"
    the same way every other soft threshold in this script does, and the
    result is combined with the percentile weight by MULTIPLICATION, so the
    transition gate can only ever ADD suppression on top of what the
    percentile test already decided, never remove it -- a pixel the
    percentile test already protects (a static subject's interior) stays
    protected regardless of what its transition count says, which matters
    because that count is a genuinely noisy statistic on its own
    (compression dithering alone gives an already-protected, perfectly
    static interior pixel a nontrivial transition count -- see
    build_background_stack()'s docstring). Two traps to watch for when
    computing this, both covered there: computing it on the SAME
    --bg-tolerance-r-filtered samples the percentile uses corrupts it in
    dense sky for the same reason --bg-tolerance-r itself needs
    --min-static-blob-area; and a pixel crossed by many DIFFERENT stars over
    a long window also racks up transitions, which is why the small-blob
    filter above runs a second time after this gate, not just once before
    it.
    --bg-transitions-lo/-hi are NOT fixed absolute rates by default -- they
    are auto-picked, EVERY RUN, as percentiles of THIS VIDEO's own measured
    open-sky transition-rate distribution, rather than a fixed number
    carried over from different footage -- see TRANSITION_LO_PERCENTILE's
    comment near the top of this file for why a fixed rate does not
    generalize across window lengths.
  * POINT-SOURCE TRANSIT GATE (--point-source-rate-lo/-hi,
    --point-source-transition-lo/-hi): the percentile test asks "is this
    pixel usually bright" -- which a FIXED BACKGROUND NEBULA (a Milky Way
    band that does not rotate with the starfield) answers exactly the same
    way a LIT, OPAQUE FOREGROUND SUBJECT does, even though they need
    opposite treatment. The two are not actually the same thing: a
    character or terrain is opaque, so a star passing "behind" it is
    simply hidden -- no point-source response ever appears at those
    pixels. A background nebula is not opaque -- stars pass IN FRONT of
    it, so the same top-hat response the detector already computes DOES
    appear there, arriving and departing as each star transits, exactly
    like it does over open sky. So the discriminator is not brightness but
    whether a point-source response transits the pixel at all:
    point_source_rate (see build_background_stack()) is the fraction of
    background-stack samples with a confirmed top-hat detection; where
    that rate sits in a plausible "occasional visitor" range, suppression
    is released even if bg_luma says "usually bright". A single static
    highlight (a fixed specular glint) would also cross thresh_hi
    sometimes, so rate alone is not enough -- point_source_transition_rate
    (how often that verdict flips) must ALSO stay low, i.e. look like a
    clean arrival-and-departure rather than a flicker; an idly-animated
    highlight whose fine detail jitters in place trips this second axis
    even when its rate alone would pass, the same way a swaying hair edge
    trips the (luma) transition gate above. The two axes are multiplied,
    and the result is combined with the percentile weight by MAXIMUM (not
    multiplication like the transition gate): this axis's whole purpose is
    to overturn a percentile "suppress" verdict, so it must be able to
    raise weight back up, not just lower it further. It runs BEFORE the
    small-blob filter and the (luma) transition gate, so a released nebula
    band is treated as ordinary sky by both of those from then on.
    Because the final composite (see main()) also lightens trails onto a
    bright background in proportion to this same static_weight, releasing
    it here is what actually lets a trail read as visible over the nebula
    band rather than just technically un-suppressed but invisibly dim
    against it -- see the comment at the final `max`/lighten blend.
  * --edge-protect-r spatially grows ("dilates") whatever suppression
    survives the two gates above a few pixels past its own photometric
    edge, because a soft anti-aliased silhouette boundary can itself read
    as "usually dark and settled" to both tests above even though it
    visually belongs to the foreground shape -- without the padding,
    warm-colored trails can still land right on the boundary and read as an
    implausible rim-light tint. This runs LAST, after the transition gate,
    not before it: dilating first would spread the transition gate's own
    deliberately tight, silhouette-hugging suppression further than
    necessary.
  * None of the above ever looks at WHERE a pixel is on screen. No masks, no
    regions, no hand-drawn boundary around "the character". This is what
    lets the same detector run on someone else's differently-composed
    footage unchanged.

SUB-PIXEL SAMPLING (--sample-dt, --measure-drift)
----------------------------------------------------
Trails are built by overlapping successive samples' rendered star
footprints. If a star moves MORE than about its own on-screen width between
two samples, the new footprint no longer overlaps the old one, and
consecutive samples read as a dashed line instead of a continuous streak.
So --sample-dt has to keep drift well UNDER a star's own width, not merely
under it -- this script targets ~TARGET_PX_PER_SAMPLE px of drift per
sample. Run --measure-drift to get this footage's own drift rate (px/s)
and a suggested --sample-dt; do not assume the default carries over to
different footage without checking.

DECAY IS SPECIFIED AS A REAL-TIME HALF-LIFE (--decay-half-life)
---------------------------------------------------------------------
The accumulator only ever multiplies by a fixed PER-SAMPLE factor each step
(`acc = max(acc * decay, layer)`), but this script asks for
--decay-half-life in SECONDS of real footage time instead, and converts
internally (per_sample_decay = 0.5 ** (sample_dt / half_life)). Reasoning:
per-sample decay and sample_dt are coupled -- the same per-sample factor
applied every 0.5s fades twice as fast, in real time, as applied every 1.0s.
Since --sample-dt gets re-tuned per footage (see above) for reasons that
have nothing to do with how long a trail should look, a raw per-sample decay
value would silently change the tail's real-world length every time
--sample-dt changed. Specifying seconds keeps the visual tail length
constant regardless of sampling rate. The script prints the resulting
per-sample factor and a couple of reference points (time to fade to
50%/10%) at startup.

SPATIAL SMOOTHING ALONG THE TRAIL DIRECTION (--trail-along-radius, --trail-along-mode, --trail-along-direction)
------------------------------------------------------------------------------------------------------------------
THIS IS THE SHIPPED FIX FOR THE DASHING PROBLEM, ENABLED BY DEFAULT
(--trail-along-radius 10, --trail-along-mode max, --trail-along-direction
field). `acc = max(acc*decay, layer)` is already, by construction, a
near-lossless max over every sample a single transit spans, so nothing
confined to one pixel's own brief transit has any further redundancy left
to exploit -- a pixel's rendered brightness is, to first order, just "the
twinkle's value at the one instant it transited this pixel, decayed", with
no further within-transit signal to recover. The redundancy that IS left
sits on a different axis: neighbouring pixels ALONG a trail are independent
draws of the same twinkle process at DIFFERENT times, so smoothing spatially
along the trail's own direction can recover a steadier brightness that no
per-pixel temporal filter can reach.

WHY DIRECTION HAS TO VARY ACROSS THE FRAME: the starfield is a rotating sky
dome, not a plane sliding sideways, so the true local drift direction is NOT
the same everywhere -- it visibly curves across the frame (compare a wide
render to see the fan-shaped trail pattern). A single global direction
(--measure-drift, run once) is cheap but is wrong by an increasing amount
toward the edges. --trail-along-direction=field (the default) instead
estimates a coarse grid of independently-estimated LOCAL directions and
bilinearly upsamples it to a per-pixel field -- see
estimate_local_drift_field()'s docstring for why each cell uses the same
whitened-FFT phase correlation --measure-drift uses (a plain spatial
cross-correlation was tried first and measurably locked onto the wrong peak
on this footage's dense star field, mistaking a neighbouring star's
near-identical footprint for the true match).

MAX, not average: an along-trail MAX only ever raises a dim sample toward a
brighter one drawn from a DIFFERENT transit -- unlike a temporal window
confined to one pixel's own transit, this genuinely adds information rather
than being redundant with the accumulator's own peak-hold. AVERAGE dilutes
a bright sample toward its dimmer neighbours instead of only ever filling a
dip toward a nearby peak, and does not help.

'field' vs 'global': a naive brightness-variance metric can PREFER 'global'
even though it looks worse -- 'global' can visibly smear several separate,
individually-curved trails near a frame corner (where its direction error
is largest) into wide, straight, overlapping streaks, which happens to read
as "less dashed" to a metric that only measures variance along a path, not
whether that path still belongs to one star. 'field' preserves each trail's
own curvature and keeps neighbouring trails visually distinct, so it is the
default despite scoring worse on that metric; 'global' is kept available
(cheaper: one drift measurement instead of a grid) for footage whose
rotation centre sits far enough outside frame that direction truly doesn't
vary much, but should be checked visually before trusting it.

MASKING, TO PROTECT THE HEAD AND UNTOUCHED SKY: a filter window that reaches
past a trail's own footprint would, left unchecked, either falsely light up
sky the star never actually touched (a MAX pulling in nearby lit territory)
or dim the trail's own peak by averaging in nearby unlit territory (an
AVERAGE) -- most dangerous right at the trail HEAD, since sky just past it
(in the direction of the star's future, un-recorded motion) is completely
unlit. along_trail_filter() guards against this with three rules: (1) a
pixel only ever receives a NEW (filtered) value if it already had SOME
signal (was part of the trail to begin with) -- pixels that were EXACTLY
zero in the unfiltered accumulator are left EXACTLY zero; (2) every tap the
filter draws from, including the pixel's own centre sample, is weighted by
that same lit/unlit test at the TAP's own position, so an average can never
be dragged down by counting an unlit tap in its denominator (a MAX needs no
equivalent protection: it never lowers a pixel's own value, so the only way
it could misbehave is by raising an untouched pixel, which rule (1) already
forbids); (3) every tap is additionally restricted to the TAIL side only,
never the HEAD side -- see below for why.

THE HEAD-OVERSHOOT PROBLEM, AND WHY RULES (1) AND (2) ALONE DO NOT PREVENT
IT: those two rules only guarantee that a pixel with EXACTLY zero original
signal stays exactly zero -- they say nothing about a pixel that already had
SOME (however small) genuine signal, which is exactly what a trail head
looks like close up: the star's own rendered footprint has real width, so a
handful of pixels right at, and a little past, the exact head position are
legitimately nonzero in the unfiltered accumulator to begin with (the
star's own blob, not a defect). A SYMMETRIC window lets such an
already-lit, just-past-the-head pixel look FORWARD (toward the head and
beyond) and find the head's own near-undecayed brightness within reach,
MAXing itself up to match it -- turning what should have been a small, soft
falloff (matching the base frame's own star size, invisible once
composited under it) into a hard, fully bright smear extending a further
radius_px past where the star actually was. Confirmed visually at low
--trail-opacity (which lets the base frame's full-brightness star stand out
against a dimmed trail): the base star ended up sitting in the MIDDLE of
its trail rather than at its tip.

THE FIX -- ONE-SIDED (TAIL-ONLY) SAMPLING: restrict every tap to the
pixel's TAIL side only (offsets in [-radius_px, 0], never positive) instead
of [-radius_px, +radius_px]. dir_y/dir_x already points toward the HEAD
(see along_trail_filter()'s own docstring for how that sign was established
and cross-checked), so this means: a pixel can still be raised toward a
brighter sample BEHIND it (older, filling the same twinkle-driven dashing
gaps this filter exists to fix), but never toward a sample AHEAD of it
(newer, toward the head) -- exactly the direction that let the symmetric
version borrow the head's brightness for pixels that shouldn't have had it.
A one-sided window only ever draws from half as many pixels as a symmetric
one at the same radius, so it has proportionally less material to fill a
dip with; doubling the radius largely recovers the difference. The default
is left at 10px rather than raised to compensate, since the head-overshoot
fix, not maximizing dashing reduction, is this default's job -- raise
--trail-along-radius yourself if you want stronger dashing reduction and
have checked it visually on your own footage first.

RETUNING THIS FOR DIFFERENT FOOTAGE -- WHERE TO START
----------------------------------------------------------
  1. `--measure-drift T` at a timestamp with visible stars -- gives you the
     footage's drift rate and a suggested --sample-dt.
  2. Crop a single frame to 100% zoom around one isolated star and count how
     many pixels its bright core spans (ignore any soft halo/bloom) -- that
     is roughly --tophat-k; it should be a little larger than the core.
  3. Render a short window and zoom to 100% on: an empty sky patch (checking
     for noise haze -- adjust --thresh-lo/--thresh-hi if present), a busy
     region near any static foreground subject (checking for trail residue
     on it -- adjust --bg-percentile/--bg-tolerance-r/--edge-protect-r if
     present), and the trail heads against the base frame's own stars
     (checking they coincide, which they should automatically by
     construction -- see the top of this docstring).
  4. Pick --decay-half-life by rendering a couple of candidate lengths and
     judging by eye whether streaks stay visually distinct or smear
     together; then make sure --trail-duration (or the lead-in you actually
     have available) covers 3-4x whatever half-life you land on.

One axis the steps above don't cover on their own: --min-static-blob-area,
--bg-tolerance-r, --edge-protect-r, and the hardcoded 1000px "too few pixels
to calibrate against" floor in main() are all specified in PIXELS and were
tuned on 1920x1080 footage, so switching to a 4K or 720p recording changes
what those numbers mean by itself, before anything about the scene's content
is even considered -- rescale them roughly proportionally to resolution as a
starting point, the same way you would re-measure --sample-dt for different
footage.

EXAMPLE INVOCATIONS
-----------------------
  # Measure this footage's own drift rate first, to pick --sample-dt:
  uv run sky_stack_v3.py --video night.mp4 --measure-drift 250

  # Then build the composite:
  uv run sky_stack_v3.py --video night.mp4 --out sky_final.png \\
      --base-time 390 --sample-dt 0.5 --decay-half-life 54
"""
import argparse
import json
import math
import subprocess
import sys
import time

import numpy as np
from PIL import Image
from scipy.ndimage import (
    distance_transform_edt,
    grey_opening,
    label,
    map_coordinates,
    maximum_filter,
    minimum_filter,
    zoom,
)

# Target drift, in pixels, between two consecutive trail-accumulation
# samples. See "SUB-PIXEL SAMPLING" in the module docstring: trails are
# built by overlapping each sample's rendered star footprint with the next,
# so this has to sit well under a star's on-screen width, not merely under
# it. Treat this as a starting point for a similarly-sized star, not a
# universal constant.
TARGET_PX_PER_SAMPLE = 0.3

# Rule-of-thumb number of decay half-lives needed for a trail to look fully
# developed rather than visibly tapering off short of its "natural" length.
# See "THE ONE HARD LIMITATION" in the module docstring.
HALF_LIVES_FOR_FULL_TAIL = 4.0

# How finely (real seconds between samples) the background stack should be
# sampled by default -- this drives BOTH the background-percentile statistic
# and the transition-count statistic (see the TRANSITION-COUNT GATE
# paragraph in the module docstring and build_background_stack() below), but
# only the latter actually NEEDS it fine. A percentile is a stable "usually
# dark vs usually lit" answer even from a few dozen widely-spaced samples;
# a transition COUNT is not -- it needs samples close enough together to
# actually catch a pixel flipping, or a periodic sway aliases down to
# looking like a single star's one-time sweep. Bounded to [BG_SAMPLES_MIN,
# BG_SAMPLES_MAX] so a very short window doesn't round down to too few
# samples for a stable percentile, and a very long one doesn't balloon
# memory: build_background_stack() keeps two full (n, h, w) uint8 stacks
# persistently, plus transient full-stack copies from np.partition() and
# np.stack() on top -- budget well above "two copies of the stack" when
# sizing --bg-samples for your own footage. See --bg-samples' help for a
# separate, unrelated tuning trap (compute-bound vs. decode-bound) that
# hand-lowering this value below the auto pick can walk into.
TARGET_BG_SAMPLE_DT = 1.5
BG_SAMPLES_MIN = 64
BG_SAMPLES_MAX = 400

# The transition-count gate's thresholds are NOT fixed absolute rates by
# default -- they are PERCENTILES of THIS VIDEO's OWN measured open-sky
# transition-rate distribution, recomputed every run. A fixed absolute rate
# does not generalize: how many transitions an ORDINARY single star's sweep
# contributes is not a universal constant -- it depends on what fraction of
# the window a star's transit across one pixel actually occupies (transit
# duration is fixed by drift rate and star size, roughly independent of
# window length, so a SHORT window makes that transit a much bigger slice
# of the sampled history, and the star's own twinkle can register as
# several crossings during it instead of the "1-2" a long window sees). A
# percentile of the SAME video's own open-sky population sidesteps this:
# whatever a single star's sweep does to the transition rate on THIS
# footage, at THIS window length, most open-sky pixels are experiencing
# very close to that same thing, so defining "unusually flickery" relative
# to that measured population -- rather than an absolute number carried
# over from different footage -- self-calibrates to the star density,
# drift rate, and window length actually in play.
TRANSITION_SKYLIKE_WEIGHT = 0.9   # percentile-only static_weight counted as "reference sky" for calibration
TRANSITION_LO_PERCENTILE = 85.0   # of the reference sky population's own transition-rate distribution
TRANSITION_HI_PERCENTILE = 99.0


# ---------------------------------------------------------------- utilities

def smoothstep(x, lo, hi):
    """Cubic ease 0->1 ramp between lo and hi; a soft threshold used
    throughout instead of a hard cutoff, so nothing gets a visible edge."""
    t = np.clip((x - lo) / (hi - lo), 0.0, 1.0)
    return t * t * (3 - 2 * t)


def probe_video(path):
    """ffprobe the source for width/height/fps/duration. Nothing about this
    script is allowed to assume 1920x1080 or any particular frame rate --
    other people's footage will not match this one's."""
    cmd = [
        "ffprobe", "-v", "error", "-select_streams", "v:0",
        "-show_entries", "stream=width,height,r_frame_rate",
        "-show_entries", "format=duration",
        "-of", "json", path,
    ]
    out = subprocess.run(cmd, capture_output=True, check=True).stdout
    info = json.loads(out)
    stream = info["streams"][0]
    w, h = int(stream["width"]), int(stream["height"])
    num, den = stream["r_frame_rate"].split("/")
    fps = float(num) / float(den) if float(den) else float(num)
    duration = float(info["format"]["duration"])
    return w, h, fps, duration


def grab_frame(video, t, w, h):
    """Decode exactly one frame at time t, at the source's native size."""
    cmd = [
        "ffmpeg", "-hide_banner", "-loglevel", "error", "-ss", str(t),
        "-i", video, "-frames:v", "1", "-an",
        "-f", "rawvideo", "-pix_fmt", "rgb24", "-",
    ]
    out = subprocess.run(cmd, capture_output=True, check=True).stdout
    if len(out) < w * h * 3:
        raise RuntimeError(f"failed to decode a frame at t={t} (got {len(out)} bytes, wanted {w*h*3})")
    return np.frombuffer(out, dtype=np.uint8).reshape(h, w, 3).copy()


def luma_of(frame_u8):
    return (0.299 * frame_u8[..., 0] + 0.587 * frame_u8[..., 1] + 0.114 * frame_u8[..., 2]).astype(np.uint8)


def _run_ffmpeg_decoder(cmd):
    """Launch a streaming ffmpeg decode with stdout piped for raw frames.
    stderr is deliberately left to inherit this script's own, rather than
    piped and never read: an unread stderr pipe can fill and deadlock a
    decode exactly the way an unread stdout pipe can (see _stop_ffmpeg_
    decoder() below), and -loglevel error means there is normally nothing on
    it anyway -- letting a real decode error through directly is also more
    useful than capturing and silently discarding it, which is what this
    script did before."""
    dec = subprocess.Popen(cmd, stdout=subprocess.PIPE)
    if dec.stdout is None:
        raise RuntimeError("ffmpeg subprocess has no stdout pipe -- was it opened with stdout=PIPE?")
    return dec


def _stop_ffmpeg_decoder(dec):
    """Close the child's stdout before waiting on it, and make sure the
    child is actually gone even if this is reached via an exception out of
    the caller's read loop.

    Closing stdout FIRST (rather than just looping until read_count reaches
    n_samples, then calling wait()) matters: an ffmpeg build whose -t/fps
    rounding emits one more frame than the read loop asked for would
    otherwise sit blocked writing that surplus frame into a now-permanently
    -full 64KB pipe forever, since the read loop is already done and nothing
    will ever read it again -- wait() would then block forever too, with no
    output and no error. Closing our end instead makes that write fail
    (EPIPE) rather than hang, so the child actually exits and wait()
    returns. Terminating first if the child is somehow still running covers
    the same case reached via an exception instead of a clean loop exit."""
    if dec.stdout is not None:
        dec.stdout.close()
    if dec.poll() is None:
        try:
            dec.terminate()
        except OSError:
            pass
    dec.wait()


def tophat_response(luma_u8, tophat_k):
    """Grayscale opening (erode-then-dilate) reproduces anything bigger than
    the kernel almost exactly; subtracting it from the original leaves only
    small bright structures. See "HOW STAR DETECTION WORKS" in the module
    docstring for why this is the whole trick behind rejecting the Milky
    Way / clouds / a lit character, with no notion of what any of those
    actually are."""
    opened = grey_opening(luma_u8, size=(tophat_k, tophat_k))
    th = luma_u8.astype(np.int16) - opened.astype(np.int16)
    np.clip(th, 0, None, out=th)
    return th.astype(np.float32)


# ----------------------------------------------------------- background stack

def _percentile_over_axis0(stack, percentile):
    """k-th order statistic nearest the requested percentile, along axis 0,
    via np.partition -- an O(n) partial sort per pixel that only needs the
    k-th order statistic, not a full sort. Shared by the background
    percentile itself and the transition-count gate's per-pixel level
    below, since both are "a percentile of a (samples, h, w) stack"."""
    n = stack.shape[0]
    k = round((percentile / 100.0) * (n - 1))
    k = max(0, min(n - 1, k))
    return np.partition(stack, k, axis=0)[k].astype(np.float32)


def build_background_stack(video, w, h, t_start, t_end, n_samples, tolerance_r, percentile, tophat_k, thresh_hi):
    """Decode n_samples evenly spaced over [t_start, t_end] ONCE, and compute
    THREE independent per-pixel statistics from that same decode -- see the
    "STATIC SUPPRESSION", "TRANSITION-COUNT GATE", and "POINT-SOURCE
    TRANSIT GATE" paragraphs in the module docstring for why each exists and
    why none can substitute for another. Returns (bg_luma, transition_rate,
    point_source_rate, point_source_transition_rate, n_samples_used).

    STATISTIC 1 -- bg_luma: the per-pixel LOW PERCENTILE of luma, the
    "usually swept-across sky vs. usually static foreground" test. Before
    folding each sample into this one, a small spatial max-filter (radius
    tolerance_r) is applied: it absorbs sub-pixel rendering jitter (temporal
    antialiasing, dithering), and -- more importantly on footage with an
    idly-animated character -- it reframes a soft anti-aliased silhouette
    edge's question from "was THIS pixel bright on this sample" to "was it
    bright ANYWHERE in this small neighborhood on this sample", which an
    edge pixel's neighborhood almost always satisfies even though the edge
    pixel itself sits on a mixed/dark gradient most of the time. Keep
    tolerance_r small relative to how far a real star drifts across the
    WHOLE window, or it starts hiding genuine motion near bright edges too.

    STATISTIC 2 -- transition_rate: for the SAME window, the fraction of
    consecutive sample-pairs where a pixel's luma crosses its own per-pixel
    midlevel (the average of that pixel's own [percentile, 100-percentile]
    range, so a dim pixel and a bright pixel are judged on their own scale,
    not a fixed brightness cutoff). Computed on the RAW, UN-tolerance-
    filtered luma -- NOT the same tolerance_r-filtered samples bg_luma uses.
    This matters: tolerance_r's max-filter makes "was some star within
    tolerance_r px at this instant" true across most of a busy sky at almost
    every sample (the same mechanism drop_small_static_blobs()'s docstring
    describes for the percentile), which inflates the transition count in
    dense sky the same way it inflates the percentile there. Using the raw
    signal avoids inheriting that trap.

    A genuinely static, well-lit interior pixel can show a HIGHER transition
    rate than a flickering silhouette edge, purely from per-frame
    compression/dithering noise sitting close to a narrow percentile band.
    This is exactly why the transition gate is combined with the percentile
    weight by MULTIPLICATION rather than used on its own: the percentile
    test already drives such a pixel's weight to ~0, so an unreliable
    transition count there is harmless -- min(anything, ~0) stays ~0. The
    transition axis is only ever DECISIVE in the band the percentile test
    leaves ambiguous. Samples also need to be close enough together in time
    to actually catch a pixel flipping, or a periodic sway aliases down to
    looking like a single star's one-time sweep -- see TARGET_BG_SAMPLE_DT.
    The transition RATE (not a raw count) is what --bg-transitions-lo/-hi
    are expressed in, specifically so the same two thresholds stay
    meaningful regardless of how many samples n_samples ends up being for a
    given window length.

    STATISTIC 3 -- point_source_rate / point_source_transition_rate: see
    "POINT-SOURCE TRANSIT GATE" in the module docstring. bg_luma answers
    "is this pixel usually bright" and cannot tell a fixed nebula glow from
    a lit character apart -- both are "usually bright". The question that
    CAN tell them apart is "does a small bright thing come and go here",
    which is exactly the top-hat response's OWN job (see tophat_response()),
    just watched over time instead of one frame. point_source_rate is the
    fraction of samples where the RAW (un-tolerance-filtered, for the same
    reason transition_rate avoids it) top-hat response exceeds thresh_hi --
    the same "confidently a point source, not noise" bar the detector
    itself uses. point_source_transition_rate is how often that yes/no
    verdict flips between consecutive samples. Both are accumulated with
    O(1)-per-sample running counters (a hit count and a previous-sample
    verdict), not a stored (n, h, w) stack, since nothing else needs the
    per-sample top-hat response once its own sample has been folded in.

    Implementation: every sample is kept TWICE in memory as an
    (n_samples, h, w) uint8 stack -- once raw, once tolerance-filtered --
    for the duration of this function only. The point-source counters add
    only a handful of (h, w) accumulators on top, not a third full stack.
    """
    span = t_end - t_start
    n_samples = max(2, n_samples)
    sample_dt = span / n_samples
    sample_fps = 1.0 / sample_dt

    cmd = ["ffmpeg", "-hide_banner", "-loglevel", "error"]
    if t_start:
        cmd += ["-ss", str(t_start)]
    cmd += ["-i", video, "-an", "-t", str(span), "-vf", f"fps={sample_fps}",
            "-f", "rawvideo", "-pix_fmt", "rgb24", "-"]
    dec = _run_ffmpeg_decoder(cmd)
    dec_stdout = dec.stdout

    frame_bytes = w * h * 3
    filt_size = 2 * tolerance_r + 1 if tolerance_r > 0 else 0
    raw_samples = []
    filt_samples = []
    psrc_hit_count = np.zeros((h, w), dtype=np.int32)
    psrc_flip_count = np.zeros((h, w), dtype=np.int32)
    psrc_prev_hit = None
    i = 0
    try:
        while i < n_samples:
            buf = dec_stdout.read(frame_bytes)
            if len(buf) < frame_bytes:
                print(f"WARNING: background stack source ran out at sample {i}/{n_samples}", file=sys.stderr)
                break
            frame = np.frombuffer(buf, dtype=np.uint8).reshape(h, w, 3)
            raw = luma_of(frame)
            raw_samples.append(raw)
            filt_samples.append(maximum_filter(raw, size=filt_size) if filt_size > 1 else raw)
            # Point-source presence for THIS sample only -- same top-hat
            # filter and thresh_hi bar the detector itself uses, on the RAW
            # (un-tolerance-filtered) luma. Folded into running counters
            # immediately rather than stored, see this function's docstring.
            psrc_hit = tophat_response(raw, tophat_k) > thresh_hi
            psrc_hit_count += psrc_hit
            if psrc_prev_hit is not None:
                psrc_flip_count += psrc_hit != psrc_prev_hit
            psrc_prev_hit = psrc_hit
            i += 1
    finally:
        # Closes the pipe (so a surplus frame can't wedge the child) and
        # waits for it, even if the loop above raised. See
        # _stop_ffmpeg_decoder()'s docstring.
        _stop_ffmpeg_decoder(dec)
    if not raw_samples:
        raise RuntimeError("background stack: no frames could be decoded from the requested window")

    raw_stack = np.stack(raw_samples, axis=0)  # (n, h, w) uint8
    del raw_samples
    filt_stack = np.stack(filt_samples, axis=0)  # (n, h, w) uint8
    del filt_samples
    n = raw_stack.shape[0]

    bg_luma = _percentile_over_axis0(filt_stack, percentile)
    del filt_stack

    # Per-pixel level for the transition test: the midpoint of this pixel's
    # OWN [percentile, 100-percentile] range, mirrored around wherever that
    # pixel's brightness happens to sit -- not a fixed luma cutoff, since a
    # dim star and a bright one both need to register a crossing on their
    # own scale. Then count how often consecutive samples land on opposite
    # sides of it: one comparison per sample, vectorized rather than a
    # literal Python loop.
    lo_t = _percentile_over_axis0(raw_stack, percentile)
    hi_t = _percentile_over_axis0(raw_stack, 100.0 - percentile)
    level = (lo_t + hi_t) / 2.0
    # NOT raw_stack.astype(np.float32) > level[None, ...]: .astype() first
    # would materialise a full separate float32 COPY of the whole (n, h, w)
    # stack (4x its uint8 bytes) before ever comparing it, on top of the
    # uint8 original and the bool result -- a substantial addition to the
    # memory budget --bg-samples' help already documents, worse at 4K.
    # Comparing the uint8 array against the float32 level directly lets
    # numpy's ufunc do the cast through its own small internal buffers
    # instead, with no full-size intermediate array; every uint8 value is
    # exactly representable in float32, so the result is numerically
    # identical to the old expression.
    above = raw_stack > level[None, ...]
    del raw_stack
    flips = np.count_nonzero(above[1:] != above[:-1], axis=0)
    transition_rate = flips.astype(np.float32) / max(1, n - 1)

    point_source_rate = psrc_hit_count.astype(np.float32) / n
    point_source_transition_rate = psrc_flip_count.astype(np.float32) / max(1, n - 1)

    return bg_luma, transition_rate, point_source_rate, point_source_transition_rate, n


def drop_small_static_blobs(static_weight, min_area, threshold=0.5):
    """Reject small connected components of "meaningfully suppressed"
    pixels (static_weight < threshold) before they ever reach the
    edge-protect dilation. See --min-static-blob-area's help for the defect
    this exists to fix: on a busy starfield, --bg-tolerance-r's own
    max-filter (radius tolerance_r, applied per SAMPLE before the
    percentile -- see build_background_stack) makes a sky pixel read as
    "elevated" in a background-stack sample whenever ANY star anywhere in
    the frame happens to drift within tolerance_r of it at that sample's
    instant -- not necessarily the same star twice, just whichever one is
    nearby right then. On a sky with many stars sweeping through over a
    long accumulation window, that is common enough, purely by chance, to
    occasionally tip an ordinary open-sky pixel's LOW percentile above
    --bg-lo -- a handful of unrelated stars each passing near the same
    pixel at different sampled instants can look statistically identical to
    one object sitting there continuously, even though the pixel itself was
    black in nearly every sample. This is a genuine defect in the
    percentile estimate, not something the edge-protect dilation
    introduces -- but edge-protect (a spatial dilation with a SQUARE
    footprint, since it is a plain scipy size=N filter) is what turns each
    such 1-4px false positive into a conspicuous flat-black square hole
    punched in open sky, because it blindly grows every sub-threshold
    pixel by the same fixed radius regardless of how large or well-attested
    the region behind it is.

    The distinguishing signal is spatial coherence, not photometry: a real
    static foreground subject (a character, terrain, a UI element) is
    suppressed as one large, spatially coherent blob, because every pixel
    on it is genuinely lit for most of the window. A busy-sky false
    positive is the opposite -- pixels that each independently got unlucky
    with a handful of different passing stars do not form a coherent
    shape, they scatter across the frame at whatever small size chance
    produced. A single area threshold well below "any plausible foreground
    subject" and well above "the largest coincidence a busy sky produces"
    separates the two cleanly without knowing anything about where either
    sits on screen.

    Implementation: label 4-connected components of (static_weight <
    threshold), then reset every component smaller than min_area back to
    1.0 (fully unsuppressed) so it never survives to the dilation step at
    all -- not shrunk, not softened, simply not treated as static.
    threshold=0.5 (the smoothstep's own midpoint) is an internal detail of
    where a component boundary is drawn, not a separate tunable knob: it
    only decides which pixels are grouped into candidate components, the
    min_area decision is what actually accepts or rejects a component,
    including ones that are only partially below 0.5 in their core.
    """
    if min_area <= 0:
        return static_weight
    mask = static_weight < threshold
    if not mask.any():
        return static_weight
    lbl, n = label(mask)
    if n == 0:
        return static_weight
    sizes = np.bincount(lbl.ravel())
    small_labels = np.flatnonzero(sizes < min_area)
    small_labels = small_labels[small_labels != 0]  # 0 is the background label, not a component
    if small_labels.size == 0:
        return static_weight
    drop_mask = np.isin(lbl, small_labels)
    out = static_weight.copy()
    out[drop_mask] = 1.0
    return out


# ------------------------------------------------------------ star detector

def star_layer(frame_u8, tophat_k, thresh_lo, thresh_hi, static_weight):
    """Compute ONE sample's contribution to the trail accumulator, as a
    continuous (never binary) weight -- see "HOW STAR DETECTION WORKS" in
    the module docstring. Returns star_layer (float32 HxWx3): the frame's
    own RGB scaled by the combined weight, ready to feed straight into
    `acc = max(acc * decay, star_layer)`.

    static_weight (HxW, 0..1) is the background-percentile-derived
    suppression term, computed ONCE for the whole run -- it does not depend
    on this sample, only on screen position, so it is passed in rather than
    recomputed per frame.
    """
    luma = luma_of(frame_u8)
    tophat = tophat_response(luma, tophat_k)
    weight = smoothstep(tophat, thresh_lo, thresh_hi) * static_weight
    return frame_u8.astype(np.float32) * weight[..., None]


# ------------------------------------------------------- drift measurement

def measure_drift(video, w, h, t, dt, tophat_k, thresh_lo):
    """Estimate starfield drift in px/s via FFT phase correlation, WITHOUT
    picking a hand-chosen patch of sky: both frames are first reduced to
    their top-hat response (see tophat_response -- the exact same filter
    the detector uses), which is ~0 everywhere except small bright points.
    A static foreground subject is large-scale and gets suppressed by the
    same opening operation, so it does not bias the correlation toward
    zero motion the way correlating raw frames would."""
    f0 = grab_frame(video, t, w, h)
    f1 = grab_frame(video, t + dt, w, h)
    m0 = tophat_response(luma_of(f0), tophat_k)
    m1 = tophat_response(luma_of(f1), tophat_k)
    m0[m0 < thresh_lo] = 0.0
    m1[m1 < thresh_lo] = 0.0

    win = np.hanning(h)[:, None] * np.hanning(w)[None, :]
    A = np.fft.fft2(m0 * win)
    B = np.fft.fft2(m1 * win)
    R = A * np.conj(B)
    R /= (np.abs(R) + 1e-8)
    r = np.fft.fftshift(np.fft.ifft2(R).real)
    py, px = np.unravel_index(np.argmax(r), r.shape)
    dy, dx = py - h // 2, px - w // 2

    def parab(v0, v1, v2):
        d = v0 - 2 * v1 + v2
        return 0.0 if d == 0 else 0.5 * (v0 - v2) / d

    if 0 < py < h - 1:
        dy += parab(r[py - 1, px], r[py, px], r[py + 1, px])
    if 0 < px < w - 1:
        dx += parab(r[py, px - 1], r[py, px], r[py, px + 1])
    dist = (dy ** 2 + dx ** 2) ** 0.5
    return dist / dt, dy / dt, dx / dt


def _unit_vectors(y, x):
    """Normalize (y, x) vector(s) to unit length; a zero vector maps to
    itself (rather than dividing by zero) since it has no meaningful
    direction to begin with."""
    n = np.hypot(y, x)
    n = np.where(n < 1e-9, 1.0, n)
    return y / n, x / n


def estimate_local_drift_field(video, w, h, t0, dt, tophat_k, thresh_lo, grid_n, win_r):
    """A coarse grid_n x grid_n grid of independently-estimated LOCAL drift
    directions, bilinearly upsampled to a per-pixel unit vector field --
    see "SPATIAL SMOOTHING ALONG THE TRAIL DIRECTION" in the module
    docstring for why a single global direction (measure_drift(), run once)
    is not good enough: the starfield is a rotating dome, and the true local
    drift direction visibly curves across the frame.

    Each grid cell re-runs the EXACT SAME phase-correlation primitive
    measure_drift() itself uses (a whitened FFT cross-power spectrum) on a
    small LOCAL window instead of the whole frame -- deliberately NOT a
    plain normalized spatial cross-correlation (patch-vs-search-region
    template matching), which measurably locks onto the WRONG match on a
    dense star field often enough to come out close to 180 degrees off (a
    neighbouring star's near-identical footprint outscoring the true
    match). Phase correlation is far more robust to this because whitening
    the cross-power spectrum (dividing by |R|) suppresses the contribution
    of any single dominant spatial frequency -- exactly what a field of
    many near-identical, near-regularly-spaced star footprints produces --
    rather than simply keeping whichever raw correlation peak happens to be
    tallest.

    Cells with too little structure to correlate (e.g. sitting over a large
    suppressed foreground region, or an empty patch of sky between stars)
    are left invalid and filled from the nearest valid cell before
    upsampling, so the field has no holes. Returns (field_y, field_x), each
    (h, w) float32 arrays of unit vectors.
    """
    f0 = grab_frame(video, t0, w, h)
    f1 = grab_frame(video, t0 + dt, w, h)
    m0 = tophat_response(luma_of(f0), tophat_k)
    m1 = tophat_response(luma_of(f1), tophat_k)
    m0[m0 < thresh_lo] = 0.0
    m1[m1 < thresh_lo] = 0.0

    ys = np.linspace(win_r + 2, h - win_r - 3, grid_n)
    xs = np.linspace(win_r + 2, w - win_r - 3, grid_n)
    dy_grid = np.full((grid_n, grid_n), np.nan)
    dx_grid = np.full((grid_n, grid_n), np.nan)
    win = np.hanning(2 * win_r)[:, None] * np.hanning(2 * win_r)[None, :]
    for i, cy in enumerate(ys):
        for j, cx in enumerate(xs):
            cyi, cxi = round(cy), round(cx)
            p0 = m0[cyi - win_r:cyi + win_r, cxi - win_r:cxi + win_r]
            p1 = m1[cyi - win_r:cyi + win_r, cxi - win_r:cxi + win_r]
            if p0.std() < 0.3 or p1.std() < 0.3:
                continue  # too little structure in this cell (e.g. suppressed foreground) -- filled from a neighbour below
            A = np.fft.fft2(p0 * win)
            B = np.fft.fft2(p1 * win)
            R = A * np.conj(B)
            R /= (np.abs(R) + 1e-8)
            rr = np.fft.fftshift(np.fft.ifft2(R).real)
            py, px = np.unravel_index(np.argmax(rr), rr.shape)
            dy_grid[i, j] = py - win_r
            dx_grid[i, j] = px - win_r

    valid = ~np.isnan(dy_grid)
    if not valid.any():
        raise RuntimeError(
            "estimate_local_drift_field: no grid cell had enough structure to correlate -- "
            "is there any visible star motion in this trail window?"
        )
    if not valid.all():
        idx = distance_transform_edt(~valid, return_distances=False, return_indices=True)
        dy_grid = dy_grid[tuple(idx)]
        dx_grid = dx_grid[tuple(idx)]

    uy, ux = _unit_vectors(dy_grid, dx_grid)
    zoom_y, zoom_x = h / grid_n, w / grid_n
    field_y = zoom(uy, (zoom_y, zoom_x), order=1)
    field_x = zoom(ux, (zoom_y, zoom_x), order=1)
    # zoom()'s output size can be off by a pixel from float rounding; pad/crop
    # back to exactly (h, w) rather than assuming it landed exactly.
    py, px = h - field_y.shape[0], w - field_y.shape[1]
    if py or px:
        field_y = np.pad(field_y, ((0, max(0, py)), (0, max(0, px))), mode="edge")[:h, :w]
        field_x = np.pad(field_x, ((0, max(0, py)), (0, max(0, px))), mode="edge")[:h, :w]
    field_y, field_x = _unit_vectors(field_y, field_x)
    return field_y.astype(np.float32), field_x.astype(np.float32)


def along_trail_filter(layer, dir_y, dir_x, radius_px, mode, step_px=0.5):
    """1-D filter (max or average) applied ALONG the local drift direction,
    to the fully-accumulated star_trails layer -- see "SPATIAL SMOOTHING
    ALONG THE TRAIL DIRECTION" in the module docstring for why this, unlike
    a purely temporal filter, actually reduces trail modulation:
    neighbouring pixels along a trail are independent draws of the star's
    own twinkle at different times, redundancy no filter confined to one
    pixel's single transit could ever reach.

    layer: (H, W, 3) float32 accumulated trail (build_star_trails()'s
    return value, BEFORE compositing against the base frame).
    dir_y, dir_x: either scalars (a single global direction) or (H, W)
    per-pixel unit-vector fields (estimate_local_drift_field()'s output).
    MUST point toward increasing time -- i.e. from a sample's TAIL side
    (older) toward its HEAD side (newer, the direction the star is
    drifting) -- for the one-sided restriction below to land on the correct
    side. Both callers of this function (main()'s 'global' branch via
    measure_drift(), and the 'field' branch via
    estimate_local_drift_field()) get this sign for free: both measure a
    displacement from a frame at time t to a frame at time t+dt (dt > 0) via
    the SAME forward-difference phase correlation, so the returned vector
    already points from the earlier position toward the later one, i.e.
    toward the head. Getting this backwards would not fail loudly: it would
    silently smear every head FURTHER forward (feeding it the exact defect
    this function exists to remove) while eating into the tail from the
    other end, so this was verified visually (100% crops at a trail head)
    rather than trusted on the math alone.
    radius_px: length, in OUTPUT PIXELS, of the ONE-SIDED kernel -- see
    below. Deliberately a pixel count, not a time span like every other
    smoothing knob in this script, because this filter operates entirely on
    the already-rendered image plane and has nothing to do with --sample-dt.
    mode: 'max' or 'avg'.

    ONE-SIDED (TAIL-ONLY) SAMPLING -- FOR THE TRAIL HEAD: a symmetric window
    (drawing taps from BOTH sides of a pixel) is wrong specifically at a
    trail's HEAD: a star's head is its most recent recorded position, so by
    construction there is nothing real ahead of it, but a symmetric window
    doesn't know that and happily reaches forward anyway. Near the head,
    that forward reach lets a pixel that is already lit (the star's own
    rendered footprint has some width, so a few pixels right at and just
    past the exact head position are legitimately nonzero to begin with --
    see the MASKING paragraph below) get MAXED UP to the head's own
    near-undecayed brightness, turning what should be a small, soft falloff
    matching the base frame's own star size into a hard, fully-bright smear
    extending a further radius_px past where the star actually was --
    confirmed visually, the base frame's own star (which is supposed to sit
    AT THE TIP of its trail, see build_star_trails()'s docstring) ended up
    sitting in the MIDDLE of one that had grown past it. Restricting every
    tap to the TAIL side (offsets in [-radius_px, 0], never > 0) removes
    exactly the direction responsible: a pixel can still be raised toward a
    brighter sample BEHIND it (older, filling the same dashing gaps this
    filter was built to fix), but never toward a sample AHEAD of it (newer,
    toward the head) -- so the head itself, the single most head-ward lit
    pixel in any given local neighbourhood, has no tail-side sample left
    that is BRIGHTER than itself to be pulled toward (it already IS the
    local peak), and nothing forward of it can borrow its brightness
    either, since forward is now off-limits to everyone. Applies to both
    'max' and 'avg' modes for consistency, even though 'max' is the
    default.

    MASKING (see the module docstring for the full reasoning): a pixel only
    ever receives a new, filtered value if the ORIGINAL layer already had
    some signal there (channel-max > 0) -- sky the accumulator never
    touched stays untouched, so the trail can never spatially grow past its
    own footprint out of nothing. Every tap the filter draws from, including
    the centre tap, is weighted by that same lit/unlit test sampled at the
    tap's own offset position, so an average is never dragged down by
    counting an unlit tap in its denominator; a max needs no equivalent
    protection since it never lowers a pixel's own value in the first
    place. This alone does NOT prevent the head-overshoot case above -- it
    only guarantees a pixel with EXACTLY zero original signal stays exactly
    zero, and the overshooting pixels were never exactly zero to begin with
    (they were the head's own small but genuine rendered footprint); the
    one-sided restriction is what actually fixes that, by removing the
    forward-reaching taps that let those already-nonzero pixels borrow far
    more brightness than their own original value.
    """
    H, W, C = layer.shape
    yy, xx = np.mgrid[0:H, 0:W].astype(np.float32)
    lit = layer.max(axis=2) > 0.0

    n_taps = max(1, round(radius_px / step_px))
    # ONE-SIDED: only the tail side (offsets <= 0). dir_y/dir_x points
    # toward the HEAD (see the sign-convention paragraph above), so a
    # NEGATIVE offset here samples toward the TAIL (older positions) and a
    # positive one would sample toward the HEAD -- exactly the direction
    # this function must never draw from. 0 (the pixel's own original
    # value) is always included as the tail-most-recent endpoint of the
    # window.
    offsets = np.linspace(-radius_px, 0.0, n_taps + 1)

    if mode == "max":
        acc = layer.copy()
    else:
        wsum = np.zeros((H, W), dtype=np.float32)
        vsum = np.zeros((H, W, C), dtype=np.float32)

    for k in offsets:
        if k == 0:
            sy, sx = yy, xx
        else:
            sy = yy + k * dir_y
            sx = xx + k * dir_x
        wk = map_coordinates(lit.astype(np.float32), [sy, sx], order=1, mode="constant", cval=0.0)
        if mode == "max":
            samp = np.empty((H, W, C), dtype=np.float32)
            for c in range(C):
                samp[..., c] = map_coordinates(layer[..., c], [sy, sx], order=1, mode="constant", cval=0.0)
            candidate = np.where(wk[..., None] > 0.5, samp, -1.0)
            np.maximum(acc, candidate, out=acc)
        else:
            for c in range(C):
                sc = map_coordinates(layer[..., c], [sy, sx], order=1, mode="constant", cval=0.0)
                vsum[..., c] += sc * wk
            wsum += wk

    if mode == "max":
        out = acc
    else:
        denom = np.where(wsum > 1e-6, wsum, 1.0)[..., None]
        out = vsum / denom

    return np.where(lit[..., None], out, layer).astype(np.float32)


# ------------------------------------------------------------- compositing

def build_star_trails(video, w, h, trail_start, trail_end, sample_dt, decay,
                       tophat_k, thresh_lo, thresh_hi, static_weight):
    """Decode [trail_start, trail_end] at a fixed sample_dt via ffmpeg's own
    fps filter (one sequential decode pass, cheaper than seeking frame by
    frame), computing each sample's continuous star_layer (see star_layer())
    and decay-accumulating as we go, OLDEST SAMPLE FIRST:

        acc = max(acc * decay, star_layer)

    trail_end is always --base-time, but the base frame AT that exact
    timestamp is NOT one of the samples decoded here. Samples land at
    trail_start + k*sample_dt for k = 0..n_samples-1, so the LAST sample
    actually fed into the accumulator sits at trail_end - sample_dt -- one
    whole --sample-dt short of the base frame, not at it. The base frame
    itself is decoded completely separately, by grab_frame() in main(), and
    is never passed through this function at all.

    Because samples are fed in chronological order and every step both
    decays the running accumulator AND lays the new sample on top at full
    strength, that last (trail_end - sample_dt) sample still ends up the
    brightest, undecayed layer this function produces, with every earlier
    sample's footprint fading geometrically behind it. The trail head landing
    on the base frame's own star is NOT this function matching them up --
    main() does that afterward via `final = max(base, trails)`: wherever the
    base frame's own full-brightness star is brighter than whatever this
    accumulator produced at that pixel, the base star simply wins outright
    and IS the visible head. That relies on this function's freshest sample
    sitting close enough (in time, hence in on-screen position) to the base
    frame's star that the two blend into one continuous streak -- true at
    this script's recommended drift of ~TARGET_PX_PER_SAMPLE per sample, but
    not guaranteed in general: a --sample-dt too coarse for the footage's
    drift rate would leave this function's own last, undecayed sample a
    visible full step behind the base frame's star instead of seamless with
    it. See the top of the module docstring.
    """
    span = trail_end - trail_start
    sample_fps = 1.0 / sample_dt
    # round() alone can legitimately land on 0 for a very short window (0 <
    # span < sample_dt/2) -- but ffmpeg's fps filter still emits at least one
    # frame for any non-empty -t window regardless of what n_samples says, so
    # a 0 here used to mean nothing would ever read that one frame back out
    # of the pipe: ffmpeg would block writing it into a full 64KB pipe
    # forever, and dec.wait() would then block forever too, with no output
    # and no error. This function is only ever called with span > 0 (see
    # main()'s `if trail_duration <= 0.0` branch), so max(1, ...) always
    # requests at least the one frame ffmpeg is going to produce anyway.
    # See also _stop_ffmpeg_decoder() for the other half of this fix (making
    # the read loop tolerate a surplus frame instead of assuming n_samples
    # is exact).
    n_samples = max(1, round(span * sample_fps))
    print(
        f"star trails: {n_samples} samples over t=[{trail_start:.2f}, {trail_end:.2f}] "
        f"({sample_fps:.4f} samples/s, {sample_dt}s apart, last sample at "
        f"t={trail_end - sample_dt:.2f} -- one sample short of the base frame, see "
        "this function's docstring)",
        file=sys.stderr,
    )

    cmd = ["ffmpeg", "-hide_banner", "-loglevel", "error"]
    if trail_start:
        cmd += ["-ss", str(trail_start)]
    cmd += ["-i", video, "-an", "-t", str(span), "-vf", f"fps={sample_fps}",
            "-f", "rawvideo", "-pix_fmt", "rgb24", "-"]
    dec = _run_ffmpeg_decoder(cmd)
    dec_stdout = dec.stdout

    acc = np.zeros((h, w, 3), dtype=np.float32)
    frame_bytes = w * h * 3
    t0 = time.time()
    i = 0
    try:
        while i < n_samples:
            buf = dec_stdout.read(frame_bytes)
            if len(buf) < frame_bytes:
                print(f"WARNING: source ran out at sample {i}/{n_samples}", file=sys.stderr)
                break
            frame = np.frombuffer(buf, dtype=np.uint8).reshape(h, w, 3)
            layer = star_layer(frame, tophat_k, thresh_lo, thresh_hi, static_weight)
            acc *= decay
            np.maximum(acc, layer, out=acc)
            i += 1
            if i % 50 == 0:
                print(f"  {i}/{n_samples} elapsed={time.time()-t0:.1f}s", file=sys.stderr)
    finally:
        # Closes the pipe (so a surplus frame can't wedge the child) and
        # waits for it, even if the loop above raised. See
        # _stop_ffmpeg_decoder()'s docstring.
        _stop_ffmpeg_decoder(dec)
    print(f"star trails done: {i} samples in {time.time()-t0:.1f}s", file=sys.stderr)
    return acc


# ------------------------------------------------------------------- main

def main():
    ap = argparse.ArgumentParser(
        description="Composite a decay-style star-trail still, ending exactly at one chosen base frame, from a static-camera night-sky video.",
        formatter_class=argparse.RawDescriptionHelpFormatter,
        epilog=__doc__,
    )
    ap.add_argument("--video", required=True, help="path to the source video")
    ap.add_argument("--out", default="sky_final.png", help="output PNG path for the composite")

    ap.add_argument("--base-time", type=float, default=None, help="required unless --measure-drift is given. timestamp (s) of the FINAL frame -- everything in the composite (the sky, any foreground subject, and where every star trail ends) is exactly what is visible at this instant, plus trailing history accumulated from footage BEFORE it. IMPORTANT: trail length is capped by how much footage precedes this timestamp (see \"THE ONE HARD LIMITATION\" in the module docstring) -- a --base-time too close to the start of the video will produce visibly short, tapering trails no matter how --decay-half-life or --trail-duration are set, because there simply isn't enough footage before it to accumulate a full tail. The script prints the available lead-in and how many decay half-lives it covers every run; check that line before judging a short-trail result as a bug")

    ap.add_argument("--trail-duration", type=float, default=390.0, help="how much footage (s) BEFORE --base-time to accumulate star trails over; the window is always [base_time - trail_duration, base_time], never anywhere else. Automatically clamped down to however much footage actually precedes --base-time (or to 0 if --base-time is at/near the start of the video, in which case the output is just the base frame with no trails). --trail-duration and --decay-half-life are independent knobs -- changing one does not move this default. The rule for CHOOSING a value is ~4x whatever --decay-half-life you land on (see HALF_LIVES_FOR_FULL_TAIL): that is the minimum for a trail to fully develop rather than visibly taper off short; past that point more --trail-duration doesn't change the visible trail, since anything older has already decayed to near-nothing. A clip with LESS lead-in than this default assumes will have the window auto-clamped down to whatever actually precedes --base-time, and if that clamped window covers fewer than ~4 half-lives the script prints a short-lead-in warning -- lower --trail-duration yourself to whatever your own footage's lead-in and chosen --decay-half-life actually support")
    ap.add_argument("--sample-dt", type=float, default=0.5, help=f"seconds between accumulation samples. Trails are built by overlapping each sample's star footprint with the next (see \"SUB-PIXEL SAMPLING\" in the module docstring), so the step must stay well UNDER a star's own on-screen width, not merely under it -- aim for drift_px_per_s * sample_dt around {TARGET_PX_PER_SAMPLE}px or less (run --measure-drift to get a suggestion for THIS footage). Too coarse and trails visibly bead or dash; too fine costs decode/accumulate time for diminishing visual gain past that point. Always re-measure for different footage rather than trusting the default")
    ap.add_argument("--decay-half-life", type=float, default=54.0, help="real-time half-life (s) of a star's trail brightness: the trail fades to 50%% after this many SECONDS of footage time, however many samples that spans. Converted internally to a per-sample factor (0.5 ** (sample_dt / half_life)) so the tail's real-world length stays constant when --sample-dt changes for an unrelated reason -- see \"DECAY IS SPECIFIED AS A REAL-TIME HALF-LIFE\" in the module docstring. Higher = longer, softer comet tails; lower = short, crisp dashes; float('inf') never fades (full arcs, needs --trail-duration to cover the whole lead-in); 0 shows only the base frame's own stars with no trail at all. Render a couple of candidates on your own footage before trusting the default blindly")
    ap.add_argument("--trail-opacity", type=float, default=0.35, help="multiplies the fully-accumulated star-trail layer by this factor before the final lighten against the base frame (final = max(base, trails * opacity)), so trails read as more transparent -- background sky showing through between them -- WITHOUT shortening them. Deliberately not alpha blending (base*(1-a) + trails*a): that would also darken the base frame's own sky wherever trails are dim or absent, which is not what a more-transparent trail means here; scaling only the trail layer before the max leaves the base frame untouched off-trail. 0.35 (default) dims the trail layer to just over a third strength so the background sky shows through between individual trails instead of the sky reading as one dense, nearly-opaque tangle of streaks. Pass 1.0 for full-strength, undimmed trails; values in between scale linearly")
    ap.add_argument("--trail-along-radius", type=float, default=10.0, help="length, in OUTPUT PIXELS (not seconds -- this filter runs on the already-rendered trail image and has nothing to do with --sample-dt), of a ONE-SIDED, TAIL-ONLY 1-D filter applied ALONG the local drift direction to the fully-accumulated star_trails layer, AFTER accumulation and BEFORE compositing against the base frame -- see \"SPATIAL SMOOTHING ALONG THE TRAIL DIRECTION\" in the module docstring for why this reduces trail dashing where a purely temporal filter cannot, and why it is one-sided (a symmetric window let a trail visibly overshoot past its own head, ending up with the base frame's own star sitting in the MIDDLE of its trail instead of at the tip; see along_trail_filter()'s docstring for the fix). Pass --trail-along-radius 0 to disable this entirely and reproduce the plain accumulated trail bit-for-bit, the same guarantee every other knob in this script makes for its own off state. This default was only validated visually on a small set of footage -- render a before/after comparison on your own footage rather than trusting it blindly. Costs one extra O(radius) pass of map_coordinates() calls over the whole frame, cheap relative to the accumulation pass itself")
    ap.add_argument("--trail-along-mode", choices=["max", "avg"], default="max", help="statistic the --trail-along-radius filter uses: 'max' (default) only ever raises a dim tap toward a nearby brighter one along the trail. 'avg' dilutes a bright tap toward its dimmer neighbours instead of only ever filling a dip toward a nearby peak, and was measured not to help -- kept available as a correctly-implemented alternative for footage that might differ. Both modes draw taps from the TAIL side only, never the HEAD side -- see along_trail_filter()'s docstring for why a symmetric window is wrong specifically at a trail's head")
    ap.add_argument("--trail-along-direction", choices=["field", "global"], default="field", help="how --trail-along-radius picks the local drift direction to steer along. 'field' (default) estimates a coarse --trail-along-grid x --trail-along-grid grid of INDEPENDENT local directions (phase correlation, same primitive as --measure-drift, just run on small local windows) and bilinearly upsamples it to a per-pixel field -- necessary because the starfield is a rotating dome, not a sliding plane, so the true drift direction visibly curves across the frame. 'global' uses one direction for the whole frame (this script's own --measure-drift, evaluated once) -- cheaper, but was found to visibly smear several separate, individually-curved trails near a frame corner into wide, straight, overlapping streaks, despite scoring better on a naive brightness-variance metric that doesn't notice the smearing. 'field' preserves each trail's own curvature and keeps neighbouring trails visually distinct. Only use 'global' after checking it visually on your own footage; 'field' is the safer default")
    ap.add_argument("--trail-along-grid", type=int, default=6, help="grid resolution (--trail-along-grid x --trail-along-grid cells) for --trail-along-direction=field's local drift-direction estimate. Finer grids track direction curvature more precisely but each cell needs enough local star density to correlate against (cells with too little structure are filled from their nearest valid neighbour, see estimate_local_drift_field()'s docstring) -- raising this on sparse footage risks more cells falling back to a neighbour's estimate rather than their own")

    g = ap.add_argument_group("star detector (see \"HOW STAR DETECTION WORKS\" in module docstring)")
    g.add_argument("--tophat-k", type=int, default=7, help="white top-hat opening kernel size (px, odd); should be a bit larger than your stars' pixel size. To measure that size: crop a single frame to 100%% zoom around one isolated star and count how many pixels its bright core spans (ignore any dim halo/bloom); that count is also a reasonable estimate of the star footprint --sample-dt's guidance is targeting")
    g.add_argument("--thresh-lo", type=float, default=3.0, help="top-hat response (0-255 luma units) where the soft noise floor starts letting a sample's response through. Not an accept/reject threshold -- a gentle ramp so sensor/compression noise near 0 doesn't slowly accumulate into a grey haze, while a real (possibly twinkling, possibly faint) star's response still passes through continuously")
    g.add_argument("--thresh-hi", type=float, default=8.0, help="top-hat response where the soft noise floor reaches full weight (1.0); above this, a sample's response passes through unchanged")

    m = ap.add_argument_group("static suppression (see \"HOW STAR DETECTION WORKS\" in module docstring)")
    m.add_argument("--bg-samples", type=int, default=None, help=f"number of evenly-spaced samples across the trail window used to build BOTH the background-percentile stack and the transition-count stack (see --bg-transitions-hi). Default: auto-picked from the trail window length, aiming for ~{TARGET_BG_SAMPLE_DT}s between samples and clamped to [{BG_SAMPLES_MIN},{BG_SAMPLES_MAX}] -- printed at startup. A percentile alone would tolerate far fewer, much-more-widely-spaced samples; the transition count would not, since it needs samples close enough together to actually catch a pixel flipping (see build_background_stack()'s docstring). Set this explicitly to override the auto choice, e.g. to bound memory on a very long trail window -- budget generously, since the persistent stacks this keeps are only part of the memory bill (np.partition() and np.stack() add their own transient full-stack copies on top). This pass is normally compute-bound, with ffmpeg decode hidden behind numpy work; setting this manually low on a wide trail window, or hitting BG_SAMPLES_MAX on a very long --trail-duration, can flip it decode-bound instead -- the auto-picker avoids this")
    m.add_argument("--bg-percentile", type=float, default=15.0, help="which percentile (0-100) of each pixel's sampled luma is used as the background/motion statistic, in place of a plain minimum. 0 reproduces a plain minimum exactly (fragile -- a single dark sample from a swaying hair strand or ear during idle animation is enough to misclassify a static pixel as swept-across sky, leaving trail residue on that detail). Too high and genuinely dark sky between infrequent, sparse stars starts reading as 'usually lit' and gets suppressed, quietly losing real trail material -- check sparse sky regions, not just busy ones, if you raise this. Also sets the [percentile, 100-percentile] range the transition-count gate mirrors its own per-pixel level around, so raising it narrows that range on both ends")
    m.add_argument("--bg-tolerance-r", type=int, default=8, help="radius (px) of a small spatial max-filter applied to each background-stack sample before it's folded into the percentile stack. Absorbs sub-pixel rendering jitter, and -- more importantly -- keeps an idly-animated character's swaying silhouette edge (which sits on a soft anti-aliased brightness gradient for a large fraction of the window, not just one rare instant) from reading as 'usually dark' to the percentile test. Keep it small relative to how far a real star drifts across the WHOLE window, or genuine motion starts hiding near bright edges too; 0 disables it. Deliberately NOT applied before the transition-count gate below -- see build_background_stack()'s docstring for why that would corrupt it in dense sky the same way skipping --min-static-blob-area would corrupt the percentile there")
    m.add_argument("--bg-lo", type=float, default=20.0, help="background-percentile luma (0-255) below which a pixel counts as fully 'usually dark elsewhere' (definitely swept-across sky) for the suppression soft-ramp")
    m.add_argument("--bg-hi", type=float, default=75.0, help="background-percentile luma above which a pixel counts as 'usually lit' (static foreground -- fully suppressed from the trail layer). Needs enough headroom above your footage's own static highlights' natural brightness variation -- measure a few of yours before lowering this")
    m.add_argument("--bg-transitions-lo", type=float, default=None, help=f"transition RATE (0-1: the fraction of consecutive background-stack samples where a pixel's luma crosses its own per-pixel midlevel) below which a pixel is 'settled' on this second, independent axis and keeps whatever --bg-percentile already decided for it. Exists because 'is this pixel dark MOST of the time' cannot tell a star's one-time sweep apart from a swaying hair edge that alternates hair/sky over and over -- both can land at a similar percentile, but only the swaying edge keeps FLIPPING. See the TRANSITION-COUNT GATE paragraph in the module docstring. Default: auto-picked EVERY RUN as the P{TRANSITION_LO_PERCENTILE:.0f} of THIS video's own measured open-sky transition-rate distribution (printed at startup) rather than a fixed number -- see TRANSITION_LO_PERCENTILE's comment for why a fixed absolute rate does not generalize across window lengths. Set this explicitly only to override the auto choice")
    m.add_argument("--bg-transitions-hi", type=float, default=None, help=f"transition rate above which a pixel is forced back into full suppression regardless of what --bg-percentile decided -- this is what actually reaches a thin, isolated silhouette lock (e.g. a dangling ponytail tip) that --edge-protect-r cannot, because that lock has no already-protected neighbour for a dilation to grow from. Ramped from --bg-transitions-lo via the same smoothstep-and-multiply pattern as every other soft threshold in this script, and combined with the percentile weight by MULTIPLICATION -- so this can only ever ADD suppression on top of the percentile test, never release a pixel that test already protects. Default: auto-picked EVERY RUN as the P{TRANSITION_HI_PERCENTILE:.0f} of this video's own open-sky transition-rate distribution, for the same reason --bg-transitions-lo is auto-picked. Set explicitly to override")
    m.add_argument("--edge-protect-r", type=int, default=3, help="radius (px) of a spatial dilation applied to the suppression weight AFTER both the photometric bg-percentile test and the transition-count gate above, so protection extends a few pixels past a static silhouette's edge instead of stopping exactly at it. A dark-colored foreground edge (e.g. hair against night sky) can read as ordinary 'usually dark, usually settled' sky to both purely-statistical tests above even though it visually belongs to the silhouette, letting star trails land right on the boundary -- especially conspicuous with warm-colored stars, which tint the edge implausibly. Too large and it eats genuine trail material in open sky next to ANY suppressed structure ('dead halo' with no trails right at the edge); 0 disables it. With the transition-count gate now doing most of the work on ANIMATED edges specifically, this can often be lowered versus a version of this script without that gate -- it still matters for edges that are static but merely anti-aliased (soft, not flickering)")
    m.add_argument("--min-static-blob-area", type=int, default=8000, help="minimum size (px) of a connected component of suppressed pixels (bg-percentile weight < 0.5) for it to actually count as static foreground; smaller components are reset to fully unsuppressed BEFORE --edge-protect-r's dilation runs. Exists because a busy starfield can, purely by chance, tip a handful of unrelated open-sky pixels over the suppression threshold -- see drop_small_static_blobs()'s docstring for the exact mechanism -- and --edge-protect-r's square dilation then blows each one up into a conspicuous flat-black square hole punched in open sky. A real static foreground subject is suppressed as one large, spatially coherent blob; the busy-sky false positives this rejects are small and scattered, so a default well above that scale and well below any plausible foreground subject separates the two. Lower it toward 0 to disable (restores the old behavior of trusting every suppressed pixel regardless of shape) if your footage's real foreground detail is small (e.g. a distant prop); raise it if you have a densely star-filled sky and still see small dark holes in open sky after this fix -- check with a suppression-weight dump (see the module docstring's retuning steps) rather than guessing")
    m.add_argument("--point-source-rate-lo", type=float, default=0.01, help="fraction (0-1) of background-stack samples with a confirmed top-hat detection (see --thresh-hi) below which a pixel shows no meaningful evidence of ANY point source ever transiting it, so the POINT-SOURCE TRANSIT GATE (module docstring) leaves --bg-percentile's verdict alone. Above this, it starts treating the pixel as possibly sky a star sweeps across even if --bg-percentile called it 'usually bright' -- e.g. a fixed background nebula/Milky Way band, which is not opaque, so stars transiting it produce the same top-hat blips open sky does. Measured on this script's own test footage: ordinary open sky and a Milky Way band both sit in the low single-digit percent; a lit character's interior can too, which is why this axis alone is not decisive -- see --point-source-transition-hi")
    m.add_argument("--point-source-rate-hi", type=float, default=0.08, help="point-source detection rate (see --point-source-rate-lo) above which this axis is fully convinced a point source transits the pixel. Measured default sits just above the P99 rate this script's own open-sky/Milky-Way-band test footage produced; raise it if you see suppression being released too readily on your own footage's sky, lower it if a real nebula band still reads as fully suppressed after checking --point-source-transition-hi is not the culprit instead")
    m.add_argument("--point-source-transition-lo", type=float, default=0.02, help="fraction (0-1) of consecutive background-stack samples where the point-source detection verdict (see --point-source-rate-lo) flips, below which this SECOND axis of the POINT-SOURCE TRANSIT GATE considers the pattern a clean, occasional arrival-and-departure (what a transiting star produces) rather than a flicker. A single star's transit crosses the threshold once or twice over the whole window regardless of how often it is detected, so this stays low for real transits; a static highlight whose fine detail jitters in place from idle animation flips far more often even at a similar detection rate, which is exactly the case --point-source-rate-lo/-hi alone cannot reject -- see the POINT-SOURCE TRANSIT GATE paragraph")
    m.add_argument("--point-source-transition-hi", type=float, default=0.06, help="point-source flip rate (see --point-source-transition-lo) above which the pattern is considered too jittery to be a clean transit, and this axis's release is withheld regardless of --point-source-rate-lo/-hi -- the two axes are multiplied, so either one can veto a release, but neither can override --bg-percentile's OWN suppression on its own (this whole gate is combined with the percentile weight by MAXIMUM, so it can only release, never suppress further). Measured default sits near the P95-P99 flip rate this script's own open-sky/Milky-Way-band footage produced; this is the axis actually responsible for keeping a lit character's animated highlights suppressed even though their point-source RATE alone can resemble a real transit -- lower it if you still see trail residue on animated foreground detail, and verify with a suppression-weight dump before trusting a looser value")

    ap.add_argument("--measure-drift", type=float, default=None, metavar="T", help="measurement mode: estimate starfield drift (px/s) starting at source time T, print it, and exit without building an image. Use the result to choose --sample-dt")
    ap.add_argument("--drift-dt", type=float, default=10.0, help="time gap (s) between the two frames used by --measure-drift; longer gives a more precise rate estimate but only holds if the drift rate is steady across that gap")
    args = ap.parse_args()

    # Cheap, video-independent argument sanity checks -- fail fast, before
    # even probing the source, on inputs that would otherwise crash deep
    # inside a helper with an unhelpful raw traceback (a bare
    # ZeroDivisionError for --sample-dt, a silent half-pixel-shifted top-hat
    # response for an even --tophat-k, NaN output for a degenerate
    # threshold pair -- see smoothstep(): lo == hi is 0/0 there).
    if args.sample_dt <= 0:
        print(f"error: --sample-dt must be > 0 (got {args.sample_dt})", file=sys.stderr)
        sys.exit(1)
    if args.tophat_k <= 0 or args.tophat_k % 2 == 0:
        print(f"error: --tophat-k must be a positive odd integer (got {args.tophat_k})", file=sys.stderr)
        sys.exit(1)
    if args.thresh_lo >= args.thresh_hi:
        print(f"error: --thresh-lo ({args.thresh_lo}) must be < --thresh-hi ({args.thresh_hi})", file=sys.stderr)
        sys.exit(1)
    if args.bg_lo >= args.bg_hi:
        print(f"error: --bg-lo ({args.bg_lo}) must be < --bg-hi ({args.bg_hi})", file=sys.stderr)
        sys.exit(1)
    if args.point_source_rate_lo >= args.point_source_rate_hi:
        print(f"error: --point-source-rate-lo ({args.point_source_rate_lo}) must be < --point-source-rate-hi ({args.point_source_rate_hi})", file=sys.stderr)
        sys.exit(1)
    if args.point_source_transition_lo >= args.point_source_transition_hi:
        print(f"error: --point-source-transition-lo ({args.point_source_transition_lo}) must be < --point-source-transition-hi ({args.point_source_transition_hi})", file=sys.stderr)
        sys.exit(1)
    if args.trail_opacity < 0:
        print(f"error: --trail-opacity must be >= 0 (got {args.trail_opacity})", file=sys.stderr)
        sys.exit(1)
    if args.trail_along_radius < 0:
        print(f"error: --trail-along-radius must be >= 0 (got {args.trail_along_radius})", file=sys.stderr)
        sys.exit(1)
    if args.trail_along_grid < 2:
        print(f"error: --trail-along-grid must be >= 2 (got {args.trail_along_grid})", file=sys.stderr)
        sys.exit(1)

    w, h, fps, duration = probe_video(args.video)
    print(f"source: {args.video} -> {w}x{h} @ {fps:.3f}fps, {duration:.2f}s", file=sys.stderr)

    # The container's reported duration is not a reliable bound on what
    # ffmpeg can actually decode at -ss: it commonly exceeds the last full
    # frame's own timestamp by roughly one frame's duration. A --base-time
    # (or --measure-drift point) right at or near `duration` used to pass
    # the old `<= duration` check here and then fail deep inside
    # grab_frame() with a raw, unhelpful RuntimeError traceback -- reproduced
    # with --base-time equal to the source's exact reported duration. Back
    # the bound off by one frame so it matches what can actually be decoded.
    frame_dt = 1.0 / fps if fps > 0 else 0.0
    max_decodable_time = max(0.0, duration - frame_dt)

    if args.measure_drift is not None:
        drift_end = args.measure_drift + args.drift_dt
        if not (0.0 <= args.measure_drift <= max_decodable_time) or not (0.0 <= drift_end <= max_decodable_time):
            print(
                f"error: --measure-drift={args.measure_drift} with --drift-dt={args.drift_dt} needs both "
                f"t={args.measure_drift} and t={drift_end:.2f} to be decodable, but the source's decodable "
                f"range is [0, {max_decodable_time:.2f}] (container duration is {duration:.2f}s, source fps "
                f"is {fps:.3f})",
                file=sys.stderr,
            )
            sys.exit(1)
        rate, dy, dx = measure_drift(args.video, w, h, args.measure_drift, args.drift_dt, args.tophat_k, args.thresh_lo)
        print(f"drift over {args.drift_dt}s starting at t={args.measure_drift}: dy={dy*args.drift_dt:.2f}px dx={dx*args.drift_dt:.2f}px")
        print(f"estimated drift rate: {rate:.4f} px/s")
        suggested_dt = TARGET_PX_PER_SAMPLE / rate if rate > 0 else float("inf")
        print(f"suggested --sample-dt for <={TARGET_PX_PER_SAMPLE}px/sample: {suggested_dt:.2f}s")
        print(
            "(this assumes your stars have a pixel footprint similar to the ~1px stars "
            "this default was tuned on -- eyeball your own footage at 100% zoom and widen "
            "the target if your stars render visibly larger)"
        )
        return

    if args.base_time is None:
        print("error: --base-time is required (unless using --measure-drift)", file=sys.stderr)
        sys.exit(1)

    if not (0.0 <= args.base_time <= max_decodable_time):
        print(
            f"error: --base-time={args.base_time} is outside the source's decodable range "
            f"[0, {max_decodable_time:.2f}] (container duration is {duration:.2f}s, but the last full frame "
            f"typically lands about one frame duration, {frame_dt:.3f}s, earlier than that at this "
            f"source's {fps:.3f}fps)",
            file=sys.stderr,
        )
        sys.exit(1)

    # --decay-half-life is specified in real SECONDS, not per-sample, so the
    # tail's real-world length stays fixed when --sample-dt changes for an
    # unrelated reason -- see "DECAY IS SPECIFIED AS A REAL-TIME HALF-LIFE"
    # in the module docstring. Converted here, once, into the per-sample
    # factor the accumulator actually multiplies by every step.
    if args.decay_half_life <= 0:
        per_sample_decay = 0.0
    elif math.isinf(args.decay_half_life):
        per_sample_decay = 1.0
    else:
        per_sample_decay = 0.5 ** (args.sample_dt / args.decay_half_life)

    # THE ONE HARD LIMITATION, made concrete: how much footage actually
    # precedes --base-time, versus how much was requested. The video is
    # assumed to start at t=0, so the available lead-in is simply
    # --base-time itself. Clamp rather than error, so a workable-but-short
    # request still renders -- but print loudly, since a truncated tail is
    # easy to mistake for a bug in the detector instead of a footage limit.
    available_lead_in = args.base_time
    trail_duration = min(args.trail_duration, available_lead_in)
    trail_start = args.base_time - trail_duration
    trail_end = args.base_time

    if args.decay_half_life > 0 and not math.isinf(args.decay_half_life):
        half_lives_available = trail_duration / args.decay_half_life
    else:
        half_lives_available = float("inf") if trail_duration > 0 else 0.0
    print(
        f"lead-in available before --base-time={args.base_time}: {available_lead_in:.1f}s "
        f"(requested --trail-duration={args.trail_duration:.1f}s, using {trail_duration:.1f}s "
        f"= {half_lives_available:.2f} decay half-lives)",
        file=sys.stderr,
    )
    if half_lives_available < HALF_LIVES_FOR_FULL_TAIL:
        print(
            f"  NOTE: fewer than {HALF_LIVES_FOR_FULL_TAIL:.0f} half-lives available -- trails will visibly "
            "taper off short of their full length rather than fading smoothly to nothing. "
            "Pick a later --base-time (more lead-in) or a shorter --decay-half-life if this "
            "isn't the look you want.",
            file=sys.stderr,
        )

    print(
        "parameters: "
        f"base_time={args.base_time} trail_window=[{trail_start:.2f},{trail_end:.2f}] "
        f"sample_dt={args.sample_dt} decay_half_life={args.decay_half_life}s "
        f"tophat_k={args.tophat_k} thresh_lo={args.thresh_lo} thresh_hi={args.thresh_hi} "
        f"bg_percentile={args.bg_percentile} bg_lo={args.bg_lo} bg_hi={args.bg_hi} "
        f"bg_transitions_lo={'auto' if args.bg_transitions_lo is None else args.bg_transitions_lo} "
        f"bg_transitions_hi={'auto' if args.bg_transitions_hi is None else args.bg_transitions_hi} "
        f"min_static_blob_area={args.min_static_blob_area} edge_protect_r={args.edge_protect_r} "
        f"point_source_rate_lo={args.point_source_rate_lo} point_source_rate_hi={args.point_source_rate_hi} "
        f"point_source_transition_lo={args.point_source_transition_lo} point_source_transition_hi={args.point_source_transition_hi} "
        f"trail_along_radius={args.trail_along_radius}px trail_along_mode={args.trail_along_mode} "
        f"trail_along_direction={args.trail_along_direction} trail_along_grid={args.trail_along_grid} "
        f"trail_opacity={args.trail_opacity}",
        file=sys.stderr,
    )
    # Reference points so it's obvious at a glance what the half-life setting
    # actually buys you in real time, without doing the math by hand.
    if 0.0 < per_sample_decay < 1.0:
        t50 = args.decay_half_life
        t10 = args.decay_half_life * math.log(0.1) / math.log(0.5)
        print(
            f"decay: per-sample factor={per_sample_decay:.5f} at sample_dt={args.sample_dt}s "
            f"-> trail fades to 50% after {t50:.1f}s, to 10% after {t10:.1f}s of real footage time",
            file=sys.stderr,
        )
    elif per_sample_decay >= 1.0:
        print("decay: half-life is infinite -- trails never fade, full arcs across the whole trail window", file=sys.stderr)
    else:
        print("decay: per-sample factor is 0 -- no trails, only the base frame's own stars", file=sys.stderr)

    print(f"grabbing base frame @ t={args.base_time}", file=sys.stderr)
    base_u8 = grab_frame(args.video, args.base_time, w, h)
    base_luma = luma_of(base_u8).astype(np.float32)
    base = base_u8.astype(np.float32)

    if trail_duration <= 0.0:
        print(
            "no lead-in footage available before --base-time -- output is just the base "
            "frame, with no trails at all",
            file=sys.stderr,
        )
        final = base
    else:
        # --bg-samples drives dt for BOTH statistics build_background_stack
        # computes; auto-pick it from the window length unless the user
        # overrode it -- see TARGET_BG_SAMPLE_DT for why a flat default
        # cannot serve both a short and a long trail window equally well.
        if args.bg_samples is None:
            bg_samples = round(trail_duration / TARGET_BG_SAMPLE_DT)
            bg_samples = max(BG_SAMPLES_MIN, min(BG_SAMPLES_MAX, bg_samples))
            print(
                f"--bg-samples not set: auto-picked {bg_samples} samples over the "
                f"{trail_duration:.1f}s trail window (~{trail_duration / bg_samples:.2f}s apart)",
                file=sys.stderr,
            )
        else:
            bg_samples = args.bg_samples

        # The background stack is built over the SAME window the trail
        # accumulator will run over -- it needs to answer "usually sky vs.
        # usually static foreground, and how often does it flip" for
        # exactly the samples the accumulator is about to see, no more and
        # no less.
        print(
            f"building background stack: {bg_samples} samples across "
            f"t=[{trail_start:.2f},{trail_end:.2f}] (percentile={args.bg_percentile}, tolerance_r={args.bg_tolerance_r})",
            file=sys.stderr,
        )
        bg_luma, transition_rate, point_source_rate, point_source_transition_rate, used = build_background_stack(
            args.video, w, h, trail_start, trail_end,
            bg_samples, args.bg_tolerance_r, args.bg_percentile,
            args.tophat_k, args.thresh_hi,
        )
        print(f"  background stack used {used} samples", file=sys.stderr)

        # Soft "this pixel is usually swept-across sky, not static
        # foreground" weight from the percentile test -- see "HOW STAR
        # DETECTION WORKS" in the module docstring.
        static_weight = 1.0 - smoothstep(bg_luma, args.bg_lo, args.bg_hi)

        # POINT-SOURCE TRANSIT GATE (module docstring): overturn the
        # percentile verdict wherever a point source is seen to transit the
        # pixel -- an occasional, clean detect-then-clear pattern, not a
        # jittery one -- since that is what a star crossing a background
        # nebula looks like, and bg_luma alone cannot tell that apart from a
        # lit, opaque foreground subject. Combined by MAXIMUM, not
        # multiplication: unlike the transition gate below, this axis's job
        # is to RAISE weight the percentile test lowered, never to lower it
        # further. Runs before the small-blob filter and the transition
        # gate so a released nebula band is ordinary sky to both from here.
        point_source_release = smoothstep(
            point_source_rate, args.point_source_rate_lo, args.point_source_rate_hi
        ) * (1.0 - smoothstep(
            point_source_transition_rate, args.point_source_transition_lo, args.point_source_transition_hi
        ))
        before_psrc = int(np.sum(static_weight < 0.5))
        static_weight = np.maximum(static_weight, point_source_release)
        after_psrc = int(np.sum(static_weight < 0.5))
        print(
            f"  point-source-transit gate: released {before_psrc - after_psrc}px of percentile suppression "
            f"where a point source appears to transit rather than sit still ({before_psrc}px -> {after_psrc}px suppressed)",
            file=sys.stderr,
        )

        # --bg-transitions-lo/-hi auto-calibration: percentiles of the
        # transition rate measured on pixels THIS SAME percentile test
        # already leans towards calling sky (weight >= TRANSITION_SKYLIKE_
        # WEIGHT) -- i.e. "what does an ordinary, non-flickering pixel's
        # transition rate look like on THIS footage, at THIS window
        # length", not a number carried over from different footage. See
        # TRANSITION_LO_PERCENTILE's comment for why a fixed absolute rate
        # does not generalize across window lengths the way this does.
        sky_like = static_weight >= TRANSITION_SKYLIKE_WEIGHT
        if sky_like.sum() < 1000:
            print(
                f"  WARNING: only {int(sky_like.sum())}px look like open sky to the percentile test -- "
                "too few to calibrate the transition gate against; falling back to the whole frame",
                file=sys.stderr,
            )
            sky_like = np.ones_like(static_weight, dtype=bool)
        reference_rate = transition_rate[sky_like]
        auto_lo = float(np.percentile(reference_rate, TRANSITION_LO_PERCENTILE))
        auto_hi = float(np.percentile(reference_rate, TRANSITION_HI_PERCENTILE))
        if auto_hi <= auto_lo:
            auto_hi = auto_lo + 1e-6
        bg_transitions_lo = args.bg_transitions_lo if args.bg_transitions_lo is not None else auto_lo
        bg_transitions_hi = args.bg_transitions_hi if args.bg_transitions_hi is not None else auto_hi
        # The auto_hi <= auto_lo guard just above only protects the
        # auto-picked pair against each other. It says nothing about a
        # user-supplied --bg-transitions-lo combined with an auto-picked
        # --bg-transitions-hi (or vice versa) -- that combination can still
        # land lo >= hi, which silently inverts the transition-gate ramp
        # (smoothstep() divides by hi - lo). Check the pair actually being
        # used, whatever its origin, before it reaches smoothstep() below.
        if bg_transitions_lo >= bg_transitions_hi:
            print(
                f"error: bg_transitions_lo ({bg_transitions_lo:.3f}) must be < bg_transitions_hi "
                f"({bg_transitions_hi:.3f}) -- got an inverted or degenerate transition-gate range. "
                "If you only set one of --bg-transitions-lo/--bg-transitions-hi explicitly, the other "
                "was auto-picked and landed on the wrong side of it; set both explicitly or neither.",
                file=sys.stderr,
            )
            sys.exit(1)
        print(
            f"  transition-rate calibration: reference sky population={int(sky_like.sum())}px, "
            f"rate median={float(np.median(reference_rate)):.3f} "
            f"P{TRANSITION_LO_PERCENTILE:.0f}={auto_lo:.3f} P{TRANSITION_HI_PERCENTILE:.0f}={auto_hi:.3f} "
            f"-> using bg_transitions_lo={bg_transitions_lo:.3f} bg_transitions_hi={bg_transitions_hi:.3f}",
            file=sys.stderr,
        )

        if args.min_static_blob_area > 0:
            before = int(np.sum(static_weight < 0.5))
            static_weight = drop_small_static_blobs(static_weight, args.min_static_blob_area)
            after = int(np.sum(static_weight < 0.5))
            print(
                f"  min-static-blob-area: dropped {before - after}px of suppression that belonged to "
                f"components smaller than {args.min_static_blob_area}px ({before}px -> {after}px suppressed)",
                file=sys.stderr,
            )

        # SECOND, INDEPENDENT axis: not "how often is this pixel dark" but
        # "how many times does it change" -- see the TRANSITION-COUNT GATE
        # paragraph in the module docstring. Combined by MULTIPLICATION, the
        # same pattern star_layer() itself uses to combine the top-hat
        # weight with static_weight, so this can only ever ADD suppression
        # on top of what the percentile test above already decided, never
        # release a pixel that test already protects.
        transition_release = 1.0 - smoothstep(transition_rate, bg_transitions_lo, bg_transitions_hi)
        before_t = int(np.sum(static_weight < 0.5))
        static_weight = static_weight * transition_release
        after_t = int(np.sum(static_weight < 0.5))
        print(
            f"  transition-gate: {after_t - before_t}px newly dropped below the suppression midpoint "
            f"(flickering rather than genuinely swept-across sky; {before_t}px -> {after_t}px suppressed)",
            file=sys.stderr,
        )

        if args.min_static_blob_area > 0:
            # Same chance-coincidence risk drop_small_static_blobs()'s
            # docstring describes for the percentile test applies here too:
            # a pixel crossed by many DIFFERENT stars over a long window
            # also racks up transitions, which can tip a handful of
            # unrelated busy-sky pixels into suppression. Run the identical
            # small-blob filter again on the COMBINED result to reject those
            # the same way -- this is why it runs twice, not once.
            before2 = int(np.sum(static_weight < 0.5))
            static_weight = drop_small_static_blobs(static_weight, args.min_static_blob_area)
            after2 = int(np.sum(static_weight < 0.5))
            print(
                f"  min-static-blob-area (post-transition-gate): dropped {before2 - after2}px "
                f"({before2}px -> {after2}px suppressed)",
                file=sys.stderr,
            )

        if args.edge_protect_r > 0:
            static_weight = minimum_filter(static_weight, size=2 * args.edge_protect_r + 1)
            print(f"  edge-protect: dilated suppression weight by {args.edge_protect_r}px", file=sys.stderr)

        star_trails = build_star_trails(
            args.video, w, h, trail_start, trail_end, args.sample_dt, per_sample_decay,
            args.tophat_k, args.thresh_lo, args.thresh_hi, static_weight,
        )

        # --trail-along-radius: spatial (not temporal) smoothing along the
        # local drift direction, applied to the fully-accumulated layer,
        # AFTER accumulation and BEFORE compositing against the base frame --
        # see "SPATIAL SMOOTHING ALONG THE TRAIL DIRECTION" in the module
        # docstring. Runs by default (default radius 10px); pass
        # --trail-along-radius 0 to skip this block entirely and reproduce
        # the old (pre-along-trail-smoothing) behaviour bit-for-bit.
        if args.trail_along_radius > 0:
            # Direction is measured once, near the middle of the trail
            # window, over as long a baseline as reasonably fits (better
            # signal-to-noise for the phase correlation) without running past
            # the decodable end of the source.
            align_dt = min(30.0, max(2.0, trail_duration / 4.0))
            align_t0 = (trail_start + trail_end) / 2.0
            if align_t0 + align_dt > max_decodable_time:
                align_t0 = max(0.0, max_decodable_time - align_dt)
            print(
                f"along-trail smoothing: radius={args.trail_along_radius}px mode={args.trail_along_mode} "
                f"direction={args.trail_along_direction} (direction measured at t={align_t0:.2f}, dt={align_dt:.2f}s)",
                file=sys.stderr,
            )
            if args.trail_along_direction == "global":
                rate, gdy, gdx = measure_drift(args.video, w, h, align_t0, align_dt, args.tophat_k, args.thresh_lo)
                dir_y, dir_x = _unit_vectors(gdy, gdx)
                print(f"  global direction: dy={dir_y:.4f} dx={dir_x:.4f} (rate={rate:.4f}px/s)", file=sys.stderr)
            else:
                # win_r (the local phase-correlation window's half-size, in
                # px) scales with how far apart the grid cells sit, so a
                # coarser or finer --trail-along-grid gets a proportionally
                # sized window instead of a size tuned to one resolution.
                # win_r is close to a full cell's own spacing (not half of
                # it), so neighbouring windows overlap by nearly their own
                # width -- trading away some locality for a much stronger
                # phase-correlation signal (more stars per window); a
                # smaller, non-overlapping win_r gives a noisier, less
                # effective direction field.
                cell_span = min(h, w) / args.trail_along_grid
                win_r = int(max(30, min(300, cell_span * 0.9)))
                dir_y, dir_x = estimate_local_drift_field(
                    args.video, w, h, align_t0, align_dt, args.tophat_k, args.thresh_lo,
                    args.trail_along_grid, win_r,
                )
                print(
                    f"  local direction field: {args.trail_along_grid}x{args.trail_along_grid} grid, "
                    f"win_r={win_r}px",
                    file=sys.stderr,
                )
            if args.trail_along_mode == "avg":
                print(
                    "  WARNING: --trail-along-mode=avg was measured to NOT reduce trail modulation -- see "
                    "\"SPATIAL SMOOTHING ALONG THE TRAIL DIRECTION\" in the module docstring. Only 'max' was "
                    "measured to help.",
                    file=sys.stderr,
                )
            star_trails = along_trail_filter(star_trails, dir_y, dir_x, args.trail_along_radius, args.trail_along_mode)

        # The whole composite: plain per-channel lightening (max), nothing
        # more. The base frame's own stars are already the trail's brightest
        # (undecayed) sample -- see build_star_trails -- so max() here is not
        # blending two different depictions of the same star, it is simply
        # picking whichever of "the base frame" or "everything leading up to
        # it" is brighter at each pixel, which for the trail itself is
        # always the trail (the base frame's sky elsewhere is what "wins"
        # over an accumulator that never lit that pixel at all).
        #
        # --trail-opacity scales star_trails HERE, before the max, rather
        # than alpha-blending base and star_trails (base*(1-a) +
        # star_trails*a). Alpha blending would darken the base frame itself
        # everywhere star_trails is dim or zero -- i.e. it would dim the
        # background sky, not just the trails, which is not what "more
        # transparent trails" is supposed to mean. Scaling only star_trails
        # and then taking max(base, ...) leaves every pixel where
        # star_trails is 0 exactly equal to base (max(base, 0) == base for
        # a non-negative base), so the background is untouched at any
        # opacity and only the trail's own brightness is dimmed.
        #
        # Pure max() has one blind spot: a plain dark night sky and a lit
        # BACKGROUND NEBULA are not the same kind of "bright" to a star
        # passing over them. Over dark sky, base is ~0, so max(base, trail)
        # already reduces to just "trail" -- the formula below changes
        # nothing there. Over a nebula/Milky Way band, base itself can be
        # brighter than a decayed, opacity-dimmed trail sample, and max()
        # then shows the nebula's OWN glow and nothing else -- a trail that
        # is there (the point-source gate above says so) but too dim to
        # ever win the max reads as no trail at all, right where the
        # module docstring's "stars pass in front of it" case lives. A
        # nebula is not opaque, so its own light and a transiting star's
        # light are two real, physically separate sources that should ADD,
        # not compete for which one "wins" -- exactly like an opaque
        # foreground subject should NOT: there, static_weight (and hence
        # star_trails itself) is already suppressed to ~0, so max+0 and
        # base+0 agree and this changes nothing for her either.
        # max(a,b) + w*min(a,b) is that: pure max at w=0 (dark sky, or any
        # pixel static_weight still calls foreground), pure addition at
        # w=1 (a confirmed, released nebula pixel), ramped by the exact
        # same static_weight the accumulation stage used -- no separate
        # signal, so nothing here can disagree with what already decided
        # star_trails' own brightness. Gated additionally by bright_bg
        # (base itself, not just its low percentile, actually bright) so a
        # compression-noise floor a unit or two above 0 cannot quietly
        # nudge ordinary dark sky away from bit-exact.
        trail_c = star_trails * args.trail_opacity
        bright_bg = smoothstep(base_luma, args.bg_lo, args.bg_hi)
        blend_w = static_weight * bright_bg
        final = np.maximum(base, trail_c) + blend_w[..., None] * np.minimum(base, trail_c)

    Image.fromarray(np.clip(final, 0, 255).astype(np.uint8)).save(args.out)
    print(f"saved {args.out} ({w}x{h})", file=sys.stderr)


if __name__ == "__main__":
    main()