Eclipse timelapse re-centring: lock the sun's disc by fitting the solar limb (works through partial phases, cloud, branches and horizon refraction)

README.md

Eclipse timelapse re-centring

Removes tracking drift from a solar-eclipse timelapse by locking the sun's disc to a fixed point in frame, working from the finished video alone — no original stills required.

Written for a DWARF mini smart-telescope capture of the 12 August 2026 eclipse (1280×720, 30 fps, 75 s, HEVC), where the mount tracked imperfectly: slow drift punctuated by discrete correction jumps, the largest a 173 px lurch.

Why not track the bright region

During an eclipse the centroid of the lit area is not the centre of the sun. As the moon advances, the centroid migrates away from the true centre — by deep partial phase it is off by a large fraction of the radius. Locking onto it bakes a slow systematic wobble into the video that follows the eclipse progression.

So this fits the solar limb instead. The outer edge of the crescent is still an arc of the sun's own circle, and the sun's angular radius is effectively constant across a session.

1. Fixed-radius gradient Hough. Every strong edge pixel votes for a centre at p + R·ĝ (the intensity gradient points inward at the solar limb). Solar-limb pixels agree on one point; lunar-limb, cloud, tree and horizon pixels scatter. The accumulator is padded so the centre can be found even when it lies outside the frame — which happens once the sun is partly below the horizon.
2. Robust refinement. Inliers are selected by radius and by radial agreement of the gradient direction, then least-squares fitted. Near the horizon, differential refraction squashes the disc vertically while leaving its horizontal extent alone, so the model is an ellipse with the horizontal semi-axis pinned at R and the vertical one free. On the source clip that ratio fell smoothly from 1.005 to 0.824 — the physically expected amount.

This is inherently immune to partial occlusion. Tree branches across the disc merely remove some limb points, and the fit uses whatever arc remains; it stayed accurate down to 36% arc coverage.

Two things that are easy to get wrong

Don't smooth the track. The mount corrects in discrete jumps. Those are real image motion, and a linear smoother smears each one across a second, leaving a visible slide. Worse, a blanket median/Hampel filter rejects them: it discarded a real 173 px lurch backed by 7,300 inlier limb points. clean.py therefore overrules a sample only when it both departs from its neighbours and rests on measurably weaker evidence than they do.

Do smooth two specific regions. Where the crescent is thin near totality, or the disc sits near the horizon, the apparent centre wanders because turbulent refraction deforms the limb. That is shape change, not rigid motion — translation cannot correct it, and following it only injects jitter.

Usage

Requires ffmpeg/ffprobe and Python with opencv-python, numpy, scipy.

export ECLIPSE_SRC=/path/to/clip.mp4

python calibrate.py "$ECLIPSE_SRC" 0 120 240   # frames where the disc is still whole
python track.py     "$ECLIPSE_SRC"             # -> track_raw.npy
python clean.py     "$ECLIPSE_SRC"             # -> track_clean.npy
python render.py    "$ECLIPSE_SRC" master      # -> eclipse_locked_master.mkv
python render.py    "$ECLIPSE_SRC" delivery    # -> eclipse_locked_720p.mp4

python verify.py      eclipse_locked_master.mkv
python phaseverify.py eclipse_locked_master.mkv

calibrate.py needs frame indices where the sun is unobstructed — check the spread it reports.

Output quality

The whole geometric change is folded into a single affine, so every output pixel is resampled exactly once. Decoding goes straight to yuv444p16le: chroma is upsampled once, in 16 bit, and RGB is never involved, so there is no colour-model round trip. Full colour range is preserved end to end (verified by an identity round trip: mean offset −0.03 units, max deviation exactly the 16→10 bit quantisation step).

Lanczos4 was chosen after measuring on the sharpest limb in the clip: it matched bicubic for sharpness and showed overshoot identical to bilinear, i.e. no ringing penalty, while being the better choice for the upscaled variant.

render.py sizes the crop from the p0.5–p99.5 shift range rather than the absolute extremes, because a few extreme frames would otherwise cost ~12% of the width across the entire clip. The overflowing frames run off the edge instead. Check before relying on this that every overflowing frame sits against a dark background — clean.py reports how many there are, and they must not fall in a segment where landscape is visible.

Results on the source clip

Measure	Before	After
Solar centre wander	±40–80 px, with jumps	median residual 0.02 px, p95 0.29 px
Frames within 1 px	—	96.9%
Frame-to-frame motion (phase correlation)	1.166 px	0.543 px

33 frames (1.5%) retain up to ~30 px of wobble, in two clusters: t=40.7–47.9 s (thin, cloud-broken crescent near maximum) and t=73.6–74.2 s (near the horizon). Median arc coverage there is 0.42 against 0.61 for the clip. Those are the deliberately smoothed stretches described above.

Limitations

• Field rotation is not corrected. An alt-az mount rotates the image over a session. Sunspot correlation suggested ~0.5° over the first 10 s, but that is at the measurement floor, and past frame 300 the correlation collapses — the spots are too faint at 3.5 Mbit/s once the moon covers the disc. Recentring does not address rotation; derotation would need feature tracking.
• Refraction flattening is measured but not undone. The disc is genuinely non-circular near the horizon. b/R is recorded in track_raw.npy if you want to restore it.
• Source compression artifacts around the limb are baked in and cannot be recovered.
• verify.py re-runs the same estimator on the output, so it cannot detect a consistently wrong fit. Use phaseverify.py alongside it, and check fits visually.

calibrate.py

"""Step 1 -- measure the solar radius from the opening (unobstructed) frames.

The sun's angular size is effectively constant across a session, so a single
radius pins the limb-fit model for the whole clip. Fitting it on frames where
the disc is still whole avoids any bias from the moon.

    python calibrate.py CLIP.mp4 [frame_indices...]
"""
import sys

import cv2
import numpy as np
from scipy.optimize import least_squares

import common
import sunfit


def radius_from_full_disc(bgr):
    """Least-squares circle fit to the limb of an unobstructed disc."""
    lum = sunfit.luminance(bgr)
    thr = 0.5 * float(lum.max())
    m = (lum > thr).astype(np.uint8)
    m = cv2.morphologyEx(m, cv2.MORPH_OPEN, np.ones((5, 5), np.uint8))
    cs, _ = cv2.findContours(m, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
    c = max(cs, key=cv2.contourArea).reshape(-1, 2).astype(np.float64)
    x, y = c[:, 0], c[:, 1]
    (cx0, cy0), r0 = cv2.minEnclosingCircle(c.astype(np.float32))
    fit = least_squares(lambda p: np.hypot(x - p[0], y - p[1]) - p[2],
                        [cx0, cy0, r0])
    return fit.x


if __name__ == "__main__":
    src = common.source()
    W, H, _, fps = common.probe(src)
    frames = [int(a) for a in sys.argv[2:] if a.isdigit()] or [0, 120, 240]

    rs = []
    for i in frames:
        cx, cy, r = radius_from_full_disc(common.grab(src, i, W, H))
        print(f"frame {i:5d}: centre=({cx:8.3f},{cy:8.3f})  R={r:7.3f}")
        rs.append(r)

    R = float(np.median(rs))
    spread = float(np.max(rs) - np.min(rs))
    print(f"\nadopted R = {R:.3f} px   (spread across frames {spread:.3f} px)")
    if spread > 1.5:
        print("WARNING: radii disagree -- are all these frames unobstructed?")
    common.save_radius(R)
    print("wrote radius.txt")

clean.py

"""Step 3 -- clean the centre track.

Two ideas do the work here.

Support-based outlier rejection, not magnitude-based. A blanket median or
Hampel filter is actively harmful on this footage: it discarded a real 173 px
mount lurch that was backed by 7,300 inlier limb points. A sample is overruled
only when it BOTH departs from its neighbours AND rests on measurably weaker
evidence than they do. Genuine motion keeps full inlier support and survives; a
mis-fit does not. Steps stay sharp, which matters because the mount corrects in
discrete jumps and a linear smoother would smear each one into a visible slide.

Selective smoothing. Where the crescent is thin or the disc sits near the
horizon, the apparent centre wanders because turbulent refraction deforms the
limb -- that is shape change, not rigid motion. Translation cannot correct it,
and following it only injects jitter, so those stretches are lightly smoothed.

    python clean.py CLIP.mp4        ->  track_clean.npy
"""
import numpy as np
from scipy.ndimage import median_filter, uniform_filter1d

import common

DEV_PX = 12.0        # departure from neighbours that triggers the support test
WEAK = 0.80          # inlier count below this fraction of neighbours = weak
SCATTER_PX = 3.0     # local scatter above which a frame counts as turbulent
ARC_MAX = 0.55       # ...and only where the limb arc is short


def fill(v, bad):
    x = v.copy()
    x[bad] = np.nan
    g = np.isfinite(x)
    return np.interp(np.arange(len(x)), np.nonzero(g)[0], x[g])


def scatter(v, win=15):
    m = median_filter(v, size=win, mode="nearest")
    return median_filter(np.abs(v - m), size=win, mode="nearest") * 1.4826


def main():
    src = common.source()
    W, H, _, _ = common.probe(src)
    a = np.load("track_raw.npy")
    _, cx, cy, b, n, rms, cover, _ = a.T
    N = len(a)

    bad = ~np.isfinite(cx) | (cover < 0.20)
    print(f"rejected for weak arc / no fit : {int(bad.sum())}")

    xf, yf, nf = fill(cx, bad), fill(cy, bad), fill(n, bad)

    mx = median_filter(xf, size=5, mode="nearest")
    my = median_filter(yf, size=5, mode="nearest")
    mn = median_filter(nf, size=5, mode="nearest")
    dev = np.hypot(xf - mx, yf - my)
    weak = nf < WEAK * mn
    out = (dev > DEV_PX) & weak
    print(f"departs from neighbours >{DEV_PX:.0f}px  : {int((dev > DEV_PX).sum())}")
    print(f"  weakly supported -> overruled : {int(out.sum())}")
    print(f"  well supported   -> kept real : {int(((dev > DEV_PX) & ~weak).sum())}")

    xc = np.where(out, mx, xf)
    yc = np.where(out, my, yf)

    turb = (np.maximum(scatter(xc), scatter(yc)) > SCATTER_PX) & (cover < ARC_MAX)
    print(f"turbulence-smoothed frames     : {int(turb.sum())} "
          f"({100 * turb.sum() / N:.1f}%)")
    xs = np.where(turb, uniform_filter1d(xc, 5, mode="nearest"), xc)
    ys = np.where(turb, uniform_filter1d(yc, 5, mode="nearest"), yc)

    np.save("track_clean.npy",
            np.vstack([a[:, 0], xs, ys, bad.astype(float)]).T)

    tx, ty = np.median(xs), np.median(ys)
    dx, dy = tx - xs, ty - ys
    print("\nshift range (px):")
    for nm, v in (("x", dx), ("y", dy)):
        print(f"  {nm}: min {v.min():+7.1f}  p0.5 {np.percentile(v, 0.5):+7.1f}"
              f"  p99.5 {np.percentile(v, 99.5):+7.1f}  max {v.max():+7.1f}")
    for lbl, q in (("full extent", 0.0), ("p0.5-p99.5", 0.5)):
        lox, hix = np.percentile(dx, q), np.percentile(dx, 100 - q)
        loy, hiy = np.percentile(dy, q), np.percentile(dy, 100 - q)
        cw = W - (max(hix, 0) - min(lox, 0))
        ch = H - (max(hiy, 0) - min(loy, 0))
        over = int(((dx < lox) | (dx > hix) | (dy < loy) | (dy > hiy)).sum())
        print(f"\n[{lbl}] valid area {cw:.0f} x {ch:.0f}   "
              f"frames needing edge fill {over} ({100 * over / N:.2f}%)")
    print("\nwrote track_clean.npy")


if __name__ == "__main__":
    main()

common.py

"""Shared configuration for the eclipse re-centring pipeline."""
import json
import os
import subprocess
import sys


def source(argv_index=1):
    """Source clip: first CLI argument, else the ECLIPSE_SRC environment variable."""
    if len(sys.argv) > argv_index and sys.argv[argv_index].lower().endswith(
            (".mp4", ".mkv", ".mov", ".m4v")):
        return sys.argv[argv_index]
    p = os.environ.get("ECLIPSE_SRC")
    if not p:
        sys.exit("pass the clip as the first argument, or set ECLIPSE_SRC")
    return p


def probe(path):
    """Returns (width, height, n_frames, fps)."""
    out = subprocess.run(
        ["ffprobe", "-v", "error", "-select_streams", "v:0", "-show_entries",
         "stream=width,height,nb_frames,r_frame_rate", "-of", "json", path],
        capture_output=True, text=True).stdout
    s = json.loads(out)["streams"][0]
    num, den = s["r_frame_rate"].split("/")
    nb = s.get("nb_frames", "N/A")
    return (int(s["width"]), int(s["height"]),
            int(nb) if str(nb).isdigit() else None, float(num) / float(den))


def frame_count(path, w, h):
    """Exact frame count, decoding if the container does not state one."""
    _, _, nb, _ = probe(path)
    if nb:
        return nb
    out = subprocess.run(
        ["ffprobe", "-v", "error", "-count_frames", "-select_streams", "v:0",
         "-show_entries", "stream=nb_read_frames", "-of", "csv=p=0", path],
        capture_output=True, text=True).stdout.strip()
    return int(out)


def save_radius(r):
    with open("radius.txt", "w") as f:
        f.write(f"{r:.6f}\n")


def radius():
    try:
        with open("radius.txt") as f:
            return float(f.read().strip())
    except FileNotFoundError:
        sys.exit("run calibrate.py first (it writes radius.txt)")


def grab(path, idx, w, h):
    """Decode a single frame by index as BGR."""
    import numpy as np
    out = subprocess.run(
        ["ffmpeg", "-v", "error", "-i", path, "-vf", f"select=eq(n\\,{idx})",
         "-vsync", "0", "-frames:v", "1", "-f", "rawvideo", "-pix_fmt", "bgr24", "-"],
        capture_output=True).stdout
    return np.frombuffer(out[: w * h * 3], np.uint8).reshape(h, w, 3).copy()

phaseverify.py

"""Step 5b -- independent stability check by frame-to-frame phase correlation.

This shares no machinery with the limb fit, so it cannot inherit a systematic
bias from it. Run it on the source and on the output and compare the medians.

Read the tails with care: phase correlation also sees the moon advancing across
the disc and, at the end, the landscape sweeping through frame while the sun is
held still. Near totality the frame is almost entirely black and the correlation
peak collapses, so pairs with a weak peak are excluded from the statistics.

    python phaseverify.py CLIP.mp4
"""
import subprocess

import cv2
import numpy as np

import common

PEAK_MIN = 0.05


def main():
    path = common.source()
    W, H, _, fps = common.probe(path)
    NF = common.frame_count(path, W, H)
    FRAME = W * H

    proc = subprocess.Popen(
        ["ffmpeg", "-v", "error", "-i", path,
         "-f", "rawvideo", "-pix_fmt", "gray", "-"],
        stdout=subprocess.PIPE, bufsize=FRAME * 8)

    win = cv2.createHanningWindow((W, H), cv2.CV_32F)
    prev, shifts = None, []
    for i in range(NF):
        buf = proc.stdout.read(FRAME)
        if len(buf) < FRAME:
            break
        g = np.frombuffer(buf, np.uint8).reshape(H, W).astype(np.float32)
        if prev is not None:
            (dx, dy), resp = cv2.phaseCorrelate(prev, g, win)
            shifts.append((i, dx, dy, resp))
        prev = g
    proc.stdout.close()
    proc.wait()

    s = np.array(shifts)
    i, dx, dy, resp = s.T
    d = np.hypot(dx, dy)
    ok = resp > PEAK_MIN
    v = d[ok]
    print(f"{path}")
    print(f"  usable pairs {int(ok.sum())}/{len(s)}")
    print(f"  frame-to-frame motion: median {np.median(v):6.3f} px   "
          f"p95 {np.percentile(v, 95):6.3f}   p99 {np.percentile(v, 99):6.3f}")


if __name__ == "__main__":
    main()

render.py

"""Step 4 -- apply the lock and encode.

Quality rules enforced here:

  * decode straight to yuv444p16le, so chroma is upsampled exactly once, in 16
    bit, and RGB is never involved -- no colour-model round trip;
  * the entire geometric change (lock shift + crop origin + any scale) is folded
    into ONE affine, so every output pixel is resampled exactly once;
  * full colour range is preserved end to end (this source is yuvj420p);
  * Lanczos4 resampling -- measured on the sharpest limb in the clip it matched
    bicubic for sharpness and showed overshoot identical to bilinear, i.e. no
    ringing penalty, while being the better choice for the upscaled variant.

    python render.py CLIP.mp4 master     ->  lossless FFV1, native cropped size
    python render.py CLIP.mp4 delivery   ->  H.265 10-bit, 16:9, 1280x720
"""
import subprocess
import sys

import cv2
import numpy as np

import common

CROP_Q = 0.5         # percentile defining the crop budget (see load_track)
OUT_MASTER = "eclipse_locked_master.mkv"
OUT_DELIVERY = "eclipse_locked_720p.mp4"


def load_track(W, H, q=CROP_Q):
    """Crop budget from the p[q]..p[100-q] shift range, not the absolute extremes.

    A handful of frames carry a mount lurch far larger than the rest; sizing the
    crop for them costs ~12% of the width across the entire clip. Those frames
    instead run off the edge. Check before relying on this that every
    overflowing frame sits against a dark background -- in the source this was
    written for, none fell in the sunset segment, where a missing edge would be
    obvious against the landscape.
    """
    t = np.load("track_clean.npy")
    cx, cy = t[:, 1], t[:, 2]
    tx, ty = np.median(cx), np.median(cy)
    dx, dy = tx - cx, ty - cy
    lox, hix = np.percentile(dx, q), np.percentile(dx, 100 - q)
    loy, hiy = np.percentile(dy, q), np.percentile(dy, 100 - q)
    x0, y0 = max(hix, 0.0), max(hiy, 0.0)
    cw = int(np.floor(W - (max(hix, 0) - min(lox, 0))))
    ch = int(np.floor(H - (max(hiy, 0) - min(loy, 0))))
    return dx, dy, x0, y0, cw - (cw % 2), ch - (ch % 2)


def warp(planes, dx, dy, x0, y0, ow, oh, scale):
    """One affine per plane: lock shift, crop origin and scale in a single resample."""
    M = np.array([[scale, 0.0, scale * (dx - x0)],
                  [0.0, scale, scale * (dy - y0)]], np.float64)
    out = np.empty((3, oh, ow), np.uint16)
    for c in range(3):
        out[c] = cv2.warpAffine(planes[c], M, (ow, oh),
                                flags=cv2.INTER_LANCZOS4,
                                borderMode=cv2.BORDER_REPLICATE).astype(np.uint16)
    return out


def encoder(mode, ow, oh, fps, path):
    cmd = ["ffmpeg", "-v", "error", "-y",
           "-f", "rawvideo", "-pix_fmt", "yuv444p16le",
           "-s", f"{ow}x{oh}", "-r", f"{fps}",
           "-color_range", "pc", "-colorspace", "bt709",
           "-color_primaries", "bt709", "-color_trc", "bt709", "-i", "-"]
    if mode == "master":
        cmd += ["-c:v", "ffv1", "-level", "3", "-coder", "1", "-context", "1",
                "-g", "1", "-slices", "24", "-slicecrc", "1",
                "-pix_fmt", "yuv444p10le"]
    else:
        cmd += ["-c:v", "libx265", "-preset", "slow", "-crf", "14",
                "-pix_fmt", "yuv420p10le", "-tag:v", "hvc1", "-x265-params",
                "range=full:colorprim=bt709:transfer=bt709:colormatrix=bt709"]
    cmd += ["-color_range", "pc", "-colorspace", "bt709",
            "-color_primaries", "bt709", "-color_trc", "bt709", path]
    return subprocess.Popen(cmd, stdin=subprocess.PIPE)


def main():
    src = common.source()
    mode = next((a for a in sys.argv[1:] if a in ("master", "delivery")), "master")
    W, H, _, fps = common.probe(src)
    NF = common.frame_count(src, W, H)
    FRAME = W * H * 2 * 3

    dx, dy, x0, y0, cw, ch = load_track(W, H)
    if mode == "master":
        ow, oh, scale, sx0, sy0 = cw, ch, 1.0, x0, y0
        out = OUT_MASTER
    else:
        w9 = min(cw, ch * 16 / 9)
        scale = 1280.0 / w9
        ow, oh = 1280, 720
        sx0 = x0 + (cw - w9) / 2.0
        sy0 = y0 + (ch - w9 * 9 / 16) / 2.0
        out = OUT_DELIVERY
    print(f"valid area {cw}x{ch} -> output {ow}x{oh} (scale {scale:.4f})")

    dec = subprocess.Popen(
        ["ffmpeg", "-v", "error", "-i", src,
         "-f", "rawvideo", "-pix_fmt", "yuv444p16le", "-"],
        stdout=subprocess.PIPE, bufsize=FRAME * 3)
    enc = encoder(mode, ow, oh, fps, out)

    for i in range(NF):
        buf = dec.stdout.read(FRAME)
        if len(buf) < FRAME:
            break
        pl = np.frombuffer(buf, "<u2").reshape(3, H, W)
        enc.stdin.write(warp(pl, dx[i], dy[i], sx0, sy0, ow, oh, scale).tobytes())
        if i % 300 == 0:
            print(f"  {i}/{NF}", flush=True)

    dec.stdout.close()
    dec.wait()
    enc.stdin.close()
    enc.wait()
    print("wrote", out)


if __name__ == "__main__":
    main()

sunfit.py

"""Solar limb fitting for eclipse footage.

The lit region's centroid is useless here: as the moon advances the centroid
migrates away from the true solar centre. Instead we fit the *solar limb* --
the outer arc of the crescent is still an arc of the sun's own circle, whose
radius is fixed for the whole clip.

Stage 1: gradient-direction Hough with the radius held fixed. Every strong
edge pixel votes for a centre at p + R*grad_hat (the intensity gradient points
inward at the solar limb). Solar-limb pixels agree on one point; lunar-limb,
cloud, horizon and tree pixels scatter.

Stage 2: pick inliers by radius and by radial agreement of the gradient, then
least-squares fit. Near the horizon differential refraction squashes the disc
vertically but leaves the horizontal extent alone, so the model is an ellipse
with the horizontal semi-axis pinned at R and the vertical one free.
"""

import numpy as np
import cv2
from scipy.optimize import least_squares


def luminance(bgr):
    """Returns a float32 brightness plane robust to the red saturation at sunset."""
    return bgr.astype(np.float32).max(axis=2)


def gradients(lum, blur=2.0):
    sm = cv2.GaussianBlur(lum, (0, 0), blur)
    gx = cv2.Sobel(sm, cv2.CV_32F, 1, 0, ksize=3)
    gy = cv2.Sobel(sm, cv2.CV_32F, 0, 1, ksize=3)
    return gx, gy, cv2.magnitude(gx, gy)


def edge_points(gx, gy, mag, frac=0.25, max_pts=60000):
    """Strong edge pixels with unit gradient vectors pointing toward brighter."""
    thr = frac * float(mag.max())
    if thr <= 0:
        return None
    ys, xs = np.nonzero(mag > thr)
    if len(xs) < 30:
        return None
    if len(xs) > max_pts:                      # keep the strongest, deterministically
        m = mag[ys, xs]
        keep = np.argpartition(-m, max_pts)[:max_pts]
        ys, xs = ys[keep], xs[keep]
    g = np.stack([gx[ys, xs], gy[ys, xs]], axis=1)
    n = np.linalg.norm(g, axis=1, keepdims=True)
    return xs.astype(np.float32), ys.astype(np.float32), g / n, mag[ys, xs]


def hough_centre(pts, R, shape, pad, hint=None, hint_r=70.0):
    """Fixed-radius gradient Hough. The centre may lie outside the frame once
    the sun is partly below the horizon, hence the padded accumulator.

    `hint` restricts the search to a window around the previous centre, which
    stops thin-crescent or cloud-broken frames from latching onto a spurious
    peak somewhere else in the accumulator.
    """
    xs, ys, u, w = pts
    h, w_img = shape
    cx, cy = xs + R * u[:, 0], ys + R * u[:, 1]
    acc, _, _ = np.histogram2d(
        cy, cx, weights=w,
        bins=[h + 2 * pad, w_img + 2 * pad],
        range=[[-pad, h + pad], [-pad, w_img + pad]],
    )
    acc = cv2.GaussianBlur(acc.astype(np.float32), (0, 0), 3.0)
    if hint is not None:
        gy_, gx_ = np.mgrid[0:acc.shape[0], 0:acc.shape[1]]
        far = (gx_ - (hint[0] + pad)) ** 2 + (gy_ - (hint[1] + pad)) ** 2 > hint_r ** 2
        acc = acc.copy()
        acc[far] = 0.0
    iy, ix = np.unravel_index(int(np.argmax(acc)), acc.shape)
    # sub-bin refinement: intensity centroid of a small window around the peak
    r = 6
    y0, y1 = max(0, iy - r), min(acc.shape[0], iy + r + 1)
    x0, x1 = max(0, ix - r), min(acc.shape[1], ix + r + 1)
    win = acc[y0:y1, x0:x1].astype(np.float64)
    win = np.clip(win - win.min(), 0, None)
    s = win.sum()
    if s > 0:
        gy_, gx_ = np.mgrid[y0:y1, x0:x1]
        iy = (win * gy_).sum() / s
        ix = (win * gx_).sum() / s
    return ix + 0.5 - pad, iy + 0.5 - pad, float(acc.max())


def refine(pts, cx, cy, R, rtol=0.06, dot_min=0.90, fit_b=True):
    """Select solar-limb inliers around a coarse centre and least-squares fit.

    Horizontal semi-axis stays pinned at R; only the centre and (optionally)
    the vertical semi-axis are free, which is the correct model for refraction.
    """
    xs, ys, u, w = pts
    dx, dy = xs - cx, ys - cy
    d = np.hypot(dx, dy)
    ok = d > 1e-6
    radial = (dx[ok] * u[ok, 0] + dy[ok] * u[ok, 1]) / d[ok]
    sel = np.zeros(len(xs), bool)
    # gradient points inward at the solar limb, so agreement is negative radial
    sel[np.nonzero(ok)[0]] = (np.abs(d[ok] - R) < rtol * R) & (radial < -dot_min)
    n = int(sel.sum())
    if n < 25:
        return None
    px, py, pw = xs[sel], ys[sel], w[sel]
    sw = np.sqrt(pw / pw.max())                # weight by edge strength

    def resid(p):
        ux, uy = (px - p[0]) / R, (py - p[1]) / (p[2] if fit_b else R)
        return sw * (np.hypot(ux, uy) - 1.0)

    p0 = [cx, cy, R]
    lo = [cx - 40, cy - 40, 0.80 * R if fit_b else R - 1e-6]
    hi = [cx + 40, cy + 40, 1.05 * R if fit_b else R + 1e-6]
    try:
        r = least_squares(resid, p0, bounds=(lo, hi), xtol=1e-8, ftol=1e-8)
    except Exception:
        return None
    rms = float(np.sqrt(np.mean((r.fun / np.maximum(sw, 1e-6)) ** 2))) * R
    # angular spread of the inlier arc: a short arc makes the fit ill-conditioned
    ang = np.arctan2(py - r.x[1], px - r.x[0])
    cover = float(np.count_nonzero(np.histogram(ang, bins=36, range=(-np.pi, np.pi))[0]) / 36.0)
    return dict(cx=float(r.x[0]), cy=float(r.x[1]), b=float(r.x[2]),
                n=n, rms=rms, cover=cover)


def fit_frame(bgr, R, pad=400, fit_b=True, blur=2.0, frac=0.25,
              hint=None, hint_r=70.0):
    lum = luminance(bgr)
    gx, gy, mag = gradients(lum, blur)
    pts = edge_points(gx, gy, mag, frac=frac)
    if pts is None:
        return None
    cx, cy, peak = hough_centre(pts, R, lum.shape, pad, hint=hint, hint_r=hint_r)
    out = refine(pts, cx, cy, R, fit_b=fit_b)
    if out is None:
        return None
    out.update(hough=(cx, cy), peak=peak)
    return out

track.py

"""Step 2 -- track the solar centre through every frame.

Two-stage search per frame. A positional prior stops ambiguous frames (thin
crescent, broken cloud) from latching onto a spurious accumulator peak, but a
tracking mount makes occasional corrections far larger than any sensible search
window -- this clip had a 173 px lurch. So a free search runs alongside the
constrained one and wins whenever it is clearly better supported.

    python track.py CLIP.mp4        ->  track_raw.npy
"""
import numpy as np

import common
import sunfit


def main():
    src = common.source()
    W, H, _, fps = common.probe(src)
    NF = common.frame_count(src, W, H)
    R = common.radius()
    FRAME = W * H * 3
    print(f"{W}x{H}, {NF} frames, {fps:.3f} fps, R={R:.3f}")

    import subprocess
    proc = subprocess.Popen(
        ["ffmpeg", "-v", "error", "-i", src,
         "-f", "rawvideo", "-pix_fmt", "bgr24", "-"],
        stdout=subprocess.PIPE, bufsize=FRAME * 4)

    rows, hint, nfree = [], None, 0
    for i in range(NF):
        buf = proc.stdout.read(FRAME)
        if len(buf) < FRAME:
            print(f"stream ended at frame {i}")
            break
        img = np.frombuffer(buf, np.uint8).reshape(H, W, 3)

        rc = sunfit.fit_frame(img, R, hint=hint, hint_r=120.0) if hint else None
        rf = sunfit.fit_frame(img, R, hint=None)
        if rc is None:
            r = rf
            nfree += rf is not None
        elif rf is not None and rf["n"] > 1.25 * rc["n"]:
            r = rf
            nfree += 1
        else:
            r = rc

        if r is None:
            rows.append((i, np.nan, np.nan, np.nan, 0, np.nan, 0.0, 0.0))
            hint = None                      # prior is stale, re-acquire globally
            continue
        hint = (r["cx"], r["cy"])
        rows.append((i, r["cx"], r["cy"], r["b"], r["n"], r["rms"],
                     r["cover"], r["peak"]))
        if i % 200 == 0:
            print(f"{i:5d}/{NF}  ({r['cx']:7.2f},{r['cy']:7.2f})  "
                  f"b/R={r['b']/R:.3f} n={r['n']:6d} arc={r['cover']:.2f}",
                  flush=True)

    proc.stdout.close()
    proc.wait()
    a = np.array(rows, float)
    np.save("track_raw.npy", a)
    good = np.isfinite(a[:, 1])
    print(f"\nframes {len(a)}   fits {int(good.sum())}   "
          f"failures {int((~good).sum())}   free-search wins {nfree}")


if __name__ == "__main__":
    main()

verify.py

"""Step 5a -- re-track the rendered output and measure residual motion.

Caveat worth keeping in mind: this re-runs the SAME estimator on the output, so
it confirms the pipeline applied what it measured, but it cannot detect a fit
that was consistently wrong -- the same bias would recur. Pair it with
phaseverify.py, which shares no machinery with the limb fit.

    python verify.py OUTPUT.mkv
"""
import subprocess
import sys

import numpy as np

import common
import sunfit


def main():
    path = common.source()
    W, H, _, _ = common.probe(path)
    NF = common.frame_count(path, W, H)
    R = common.radius()
    FRAME = W * H * 3
    print(f"{path}: {W}x{H}, {NF} frames, R={R:.2f}")

    proc = subprocess.Popen(
        ["ffmpeg", "-v", "error", "-i", path,
         "-f", "rawvideo", "-pix_fmt", "bgr24", "-"],
        stdout=subprocess.PIPE, bufsize=FRAME * 4)

    rows, hint = [], None
    for i in range(NF):
        buf = proc.stdout.read(FRAME)
        if len(buf) < FRAME:
            break
        img = np.frombuffer(buf, np.uint8).reshape(H, W, 3)
        r = sunfit.fit_frame(img, R, hint=hint, hint_r=120.0 if hint else 1e9)
        if r is None:
            rows.append((i, np.nan, np.nan, 0.0))
            continue
        hint = (r["cx"], r["cy"])
        rows.append((i, r["cx"], r["cy"], r["cover"]))
    proc.stdout.close()
    proc.wait()

    a = np.array(rows, float)
    _, cx, cy, cov = a.T
    ok = np.isfinite(cx) & (cov >= 0.35)
    mx, my = np.median(cx[ok]), np.median(cy[ok])
    dx, dy = cx[ok] - mx, cy[ok] - my
    d = np.hypot(dx, dy)
    print(f"\nreliable frames {int(ok.sum())}/{len(a)}")
    print(f"residual from median centre:")
    print(f"  radial median {np.median(d):5.2f} px   p95 {np.percentile(d, 95):5.2f}"
          f"   max {d.max():6.2f}")
    for t in (0.25, 0.5, 1, 2, 5):
        m = d > t
        print(f"  > {t:4.2f} px : {int(m.sum()):4d} frames "
              f"({100 * m.sum() / len(a):.2f}%)")
    np.save("verify.npy", a)


if __name__ == "__main__":
    main()