Sylpheed/tools/re-capture/ring_angular.py

#!/usr/bin/env python3
"""Is the focus ring ROTATING, or just pulsing in brightness?

The temporal-std map of a focused button is an annulus, which both hypotheses
predict: a travelling bright feature varies every annulus pixel, and so does a
uniform fade.  Two observables separate them, and this script reports both.

  (1) TOTAL annulus brightness per frame.  A rotation moves brightness around
      the annulus and conserves the sum; an alpha pulse does not.
  (2) The 360-bin ANGULAR PROFILE, cross-correlated between frames.  A rotation
      shifts the profile by a lag; a pulse scales it in place.

CONTROL FIRST.  The angular estimator is run over a known synthetic rotation of
the run's own first frame (30/90/180/270 deg) and must recover it; the corpus
already has a centroid estimator that fails this by up to 19.8 deg, and that is
why one is not used here.

Usage:  ring_angular.py CX CY [FRAME ...]        (CX,CY in GAME coordinates)
"""
import os, sys
import numpy as np
from PIL import Image

DY, DX = 45, 1          # game(0,0) -> grab, measured by focus_ring_report.py
R_IN, R_OUT = 8.0, 18.0  # annulus radii, in px, read off the std map
NBINS = 360


def ndrotate(img, deg):
    """Bilinear rotation about the patch centre -- the control's known-positive."""
    h, w = img.shape
    cy, cx = (h - 1) / 2.0, (w - 1) / 2.0
    yy, xx = np.mgrid[0:h, 0:w]
    t = np.radians(deg)
    ys = (yy - cy) * np.cos(t) - (xx - cx) * np.sin(t) + cy
    xs = (yy - cy) * np.sin(t) + (xx - cx) * np.cos(t) + cx
    y0 = np.floor(ys).astype(int); x0 = np.floor(xs).astype(int)
    fy = ys - y0; fx = xs - x0
    out = np.zeros_like(img)
    for dy_, dx_, wgt in ((0, 0, (1 - fy) * (1 - fx)), (0, 1, (1 - fy) * fx),
                          (1, 0, fy * (1 - fx)), (1, 1, fy * fx)):
        yi = np.clip(y0 + dy_, 0, h - 1); xi = np.clip(x0 + dx_, 0, w - 1)
        ok = (y0 + dy_ >= 0) & (y0 + dy_ < h) & (x0 + dx_ >= 0) & (x0 + dx_ < w)
        out += np.where(ok, img[yi, xi] * wgt, 0.0)
    return out


def ndrotate(img, deg):
    """Bilinear rotation about the patch centre -- the control's known-positive."""
    h, w = img.shape
    cy, cx = (h - 1) / 2.0, (w - 1) / 2.0
    yy, xx = np.mgrid[0:h, 0:w]
    t = np.radians(deg)
    ys = (yy - cy) * np.cos(t) - (xx - cx) * np.sin(t) + cy
    xs = (yy - cy) * np.sin(t) + (xx - cx) * np.cos(t) + cx
    y0 = np.floor(ys).astype(int); x0 = np.floor(xs).astype(int)
    fy = ys - y0; fx = xs - x0
    out = np.zeros_like(img)
    for dy_, dx_, wgt in ((0, 0, (1 - fy) * (1 - fx)), (0, 1, (1 - fy) * fx),
                          (1, 0, fy * (1 - fx)), (1, 1, fy * fx)):
        yi = np.clip(y0 + dy_, 0, h - 1); xi = np.clip(x0 + dx_, 0, w - 1)
        ok = (y0 + dy_ >= 0) & (y0 + dy_ < h) & (x0 + dx_ >= 0) & (x0 + dx_ < w)
        out += np.where(ok, img[yi, xi] * wgt, 0.0)
    return out


def patch(path, cx, cy, half=28):
    a = np.array(Image.open(path).convert("RGB")).astype(np.float32)
    g = 0.299 * a[..., 0] + 0.587 * a[..., 1] + 0.114 * a[..., 2]
    return g[cy + DY - half:cy + DY + half, cx + DX - half:cx + DX + half]


def polar(p):
    """(total annulus brightness, 360-bin mean profile) of one patch."""
    h, w = p.shape
    yy, xx = np.mgrid[0:h, 0:w]
    cy, cx = (h - 1) / 2.0, (w - 1) / 2.0
    r = np.hypot(yy - cy, xx - cx)
    m = (r >= R_IN) & (r <= R_OUT)
    th = (np.degrees(np.arctan2(yy - cy, xx - cx)) + 360.0) % 360.0
    idx = np.clip((th[m] / 360.0 * NBINS).astype(int), 0, NBINS - 1)
    v = p[m]
    prof = np.zeros(NBINS); cnt = np.zeros(NBINS)
    np.add.at(prof, idx, v); np.add.at(cnt, idx, 1.0)
    prof = np.where(cnt > 0, prof / np.maximum(cnt, 1), np.nan)
    prof = np.nan_to_num(prof, nan=np.nanmean(prof))
    return float(v.sum()), prof


def lag(p0, p1):
    """Circular cross-correlation lag in degrees taking p0 -> p1."""
    a = p0 - p0.mean(); b = p1 - p1.mean()
    c = np.fft.irfft(np.fft.rfft(b) * np.conj(np.fft.rfft(a)), NBINS)
    k = int(np.argmax(c))
    peak = c[k] / np.sqrt((a * a).sum() * (b * b).sum())
    return (k if k <= 180 else k - 360), float(peak)


def main():
    cx, cy = int(sys.argv[1]), int(sys.argv[2])
    frames = sys.argv[3:]
    p0 = patch(frames[0], cx, cy)

    print("=== CONTROL: recover a known synthetic rotation of frame 0 ===")
    ok = True
    for deg in (30, 90, 180, 270):
        rot = ndrotate(p0, -deg)
        _, pr = polar(rot); _, pa = polar(p0)
        d, pk = lag(pa, pr)
        err = ((d - deg + 180) % 360) - 180
        flag = "ok " if abs(err) <= 3 else "FAIL"
        if abs(err) > 3:
            ok = False
        print(f"  {flag} applied {deg:4d} deg -> recovered {d:5d} deg  "
              f"(err {err:+4d}, peak {pk:.3f})")
    # negative control: a ring-free patch of the same frame must not correlate
    off = patch(frames[0], cx + 160, cy)
    _, po = polar(off); _, pa = polar(p0)
    _, pk = lag(pa, po)
    print(f"  ring-free patch of the same frame: peak {pk:.3f} (must be low)")
    if not ok:
        print("\nCONTROL FAILED — the estimator cannot measure this; stopping.")
        return 1
    print("  CONTROL PASSED\n")

    print("=== MEASUREMENT: successive live frames of the same focused ring ===")
    print(f"{'frame':<24} {'annulus sum':>12} {'vs f0 %':>9} {'lag vs f0':>10} {'peak':>7}")
    base_s, base_p = polar(p0)
    for f in frames:
        s, pr = polar(patch(f, cx, cy))
        d, pk = lag(base_p, pr)
        print(f"{os.path.basename(f):<24} {s:12.1f} {100*s/base_s:8.1f}% "
              f"{d:9d}d {pk:7.3f}")
    return 0


if __name__ == "__main__":
    sys.exit(main())