CLAUDE: fopen_turn_finder.py — use _dt rates, add glitch filter

Finite-diff vel1 was dominated by single-frame pose glitches from
algorithm failures (low altitude / dense canopy), giving 600+ deg/s
for sequences where _dt shows ~10 deg/s.

Changes:
- parse _dt entries (algorithm's smooth polynomial derivatives) as
  primary velocity metric — unaffected by pose glitches
- keep finite-diff as secondary; glitch = fd_max / dt_max
  (~1 = clean data, >>1 = pose jumps / algorithm failure)
- --max_glitch filter hides suspect sequences from console output
  (TSV always contains all rows for Calc analysis)
- removed windowed finite-diff (vel5) — _dt makes it redundant
Co-authored-by: Claude <claude@elphel.com>

scripts/fopen_turn_finder.py

@@ -20,26 +20,40 @@ A,T,R axis convention note (Elphel / imagej-elphel):
   So to find the sharpest change of course, rank by R.
   All three axes are reported so angular dynamics are fully visible.
 
-Metrics computed per sequence:
-  range      — max minus min angle over the full sequence (deg)
-  vel1       — max |Δangle| between consecutive frames / Δt  (deg/s)
-               Fast but sensitive to single-frame glitches in the data.
-  velW       — max |Δangle| over a W-frame window / Δt  (deg/s, default W=5)
-               More robust: averages over several frames, suppresses outliers.
+Velocity metrics:
+  _dt-based (PRIMARY — reliable):
+    The XML stores per-frame *_dt entries: the algorithm's own smooth
+    polynomial derivative (dx/dt, dy/dt, dz/dt, da/dt, dt/dt, dr/dt).
+    These are unaffected by single-frame pose glitches and represent
+    the algorithm's best estimate of instantaneous velocity.
+    max_dt_r/t/a = max |angular rate| in deg/s from _dt entries.
+
+  finite-diff (SECONDARY — glitch indicator):
+    Consecutive finite differences of the pose values, normalised by Δt.
+    For clean data these agree with _dt.  Large disagreement (ratio >> 1)
+    indicates pose glitches — algorithm failed, typically at low altitude
+    with dense canopy where feature matching breaks down.
+
+  glitch_ratio = max_fd_vel / max_dt_vel  (per axis).
+    ~1   → clean data, no glitches
+    >>1  → pose jumps present; _dt is still usable but sequence is suspect
+
+  range = max − min angle over the full sequence (deg).
+    Large range confirms a genuine manoeuvre rather than vibration.
 
 Note on sequence length:
   Number of frames per sequence (~40–500) varies with flight altitude
   (~25 m to ~250 m); lower altitude → smaller overlap window → shorter
-  sequences.  Because velocities are normalised by elapsed time (Δt in
-  seconds), the metrics are comparable across all sequence lengths.
+  sequences.  All velocity metrics are Δt-normalised and are comparable
+  across sequence lengths.
 
 Usage:
   python3 scripts/fopen_turn_finder.py \\
-      --root  /home/elphel/lwir16-proc/NC/linked/linked_1763232117-1763234145-v88 \\
-      [--window 5]          # frame offset for windowed velocity (default 5)
-      [--top   20]          # how many top sequences to print per axis (default 20)
+      --root   /path/to/linked_1763232117-1763234145-v88 \\
+      [--top   20]            # top sequences to print per axis (default 20)
       [--out   turns.tsv]     # summary TSV (optional)
-      [--sort  R_vel1]      # column to sort by (default R_vel1)
+      [--sort  R_dt]          # sort column (default R_dt)
+      [--max_glitch 5.0]      # hide sequences with glitch_ratio above this
 """
 
 import argparse
@@ -50,13 +64,19 @@ import numpy as np
 import csv
 
 
-# ────────────────────────────── XML parser ───────────────────────────────── #
+# ──────────────────────────── XML parser ─────────────────────────────────── #
 
 def parse_xml(path):
     """
     Parse one *-INTERFRAME.corr-xml.
-    Returns array of shape (N, 7): [timestamp, x, y, z, a, t, r]
-    sorted by timestamp.  Returns None if fewer than 2 scenes found.
+
+    Returns (pose_arr, dt_arr) where each is shape (N, 7):
+        [timestamp, x, y, z, a, t, r]
+    sorted by timestamp.  Returns (None, None) if fewer than 2 scenes found.
+
+    pose_arr  — absolute scene positions/orientations
+    dt_arr    — per-frame velocity estimates (algorithm's smooth derivatives)
+                dx/dt dy/dt dz/dt  da/dt dt/dt dr/dt  (m/s and rad/s)
     """
     PREFIX    = "EYESIS_DCT_AUX."
     SCENE_KEY = PREFIX + "scenes_"
@@ -64,15 +84,19 @@ def parse_xml(path):
         tree    = ET.parse(path)
         entries = {e.get("key"): e.text for e in tree.getroot().findall("entry")}
     except Exception:
-        return None
+        return None, None
 
-    raw = {}
+    raw_pose = {}
+    raw_dt   = {}
     for key, value in entries.items():
         if not key.startswith(SCENE_KEY):
             continue
         ts_str = key[len(SCENE_KEY):]
-        if ts_str.endswith("_dt") or ts_str.endswith("_d2t"):
+        is_dt  = ts_str.endswith("_dt")
+        if ts_str.endswith("_d2t"):
             continue
+        if is_dt:
+            ts_str = ts_str[:-3]
         try:
             vals = [float(v) for v in value.split(",")]
         except (TypeError, ValueError):
@@ -80,35 +104,43 @@ def parse_xml(path):
         if len(vals) != 6:
             continue
         # "1763233715_106222"  →  1763233715.106222 s
+        # Replace only the first underscore; the microsecond part has no underscores.
         ts = float(ts_str.replace("_", ".", 1))
-        raw[ts] = vals
+        if is_dt:
+            raw_dt[ts]   = vals
+        else:
+            raw_pose[ts] = vals
 
-    if len(raw) < 2:
-        return None
+    if len(raw_pose) < 2:
+        return None, None
 
-    rows = sorted(raw.items())
-    arr  = np.array([[ts] + vals for ts, vals in rows])  # (N,7)
-    return arr
+    def to_arr(d):
+        rows = sorted(d.items())
+        return np.array([[ts] + vals for ts, vals in rows])
 
+    pose_arr = to_arr(raw_pose)
+    dt_arr   = to_arr(raw_dt) if len(raw_dt) >= 2 else None
+    return pose_arr, dt_arr
 
-# ──────────────────────────── metrics ────────────────────────────────────── #
 
-def sequence_metrics(arr, window):
+# ─────────────────────────── metrics ─────────────────────────────────────── #
+
+def sequence_metrics(pose_arr, dt_arr):
     """
-    Compute metrics for one sequence array (N,7): [ts,x,y,z,a,t,r].
+    Compute per-sequence metrics.
 
     Returns dict with keys:
-      n_frames, duration_s,
-      for each axis in (a, t, r):
-        {axis}_range(deg)
-        {axis}_vel1_max(deg/s)   — max |Δ| between consecutive frames / Δt
-        {axis}_vel1_ts           — timestamp of that frame
-        {axis}_velW_max(deg/s)   — max |Δ| over window frames / Δt
-        {axis}_velW_ts           — timestamp of start frame
+      n_frames, duration_s
+      For each axis in (a, t, r):
+        {ax}_range(deg)
+        {ax}_dt_max(deg/s)    — max |_dt angular rate|  (PRIMARY, reliable)
+        {ax}_dt_ts            — timestamp of that frame
+        {ax}_fd_max(deg/s)    — max consecutive finite-diff angular rate
+        {ax}_glitch           — fd_max / dt_max  (1=clean, >>1=glitches)
     """
-    ts   = arr[:, 0]
-    vals = np.degrees(arr[:, 4:7])   # a, t, r  in degrees
-    n    = len(arr)
+    ts   = pose_arr[:, 0]
+    pose = np.degrees(pose_arr[:, 4:7])   # a, t, r in degrees
+    n    = len(pose_arr)
 
     result = {
         'n_frames':   n,
@@ -116,31 +148,47 @@ def sequence_metrics(arr, window):
     }
 
     axis_names = ('a', 't', 'r')
-    for ci, axis in enumerate(axis_names):
-        col = vals[:, ci]
-
-        # range
-        result[f'{axis}_range(deg)'] = float(col.max() - col.min())
-
-        # consecutive velocity  (frame i → i+1)
-        dt1   = np.diff(ts)
-        dt1   = np.where(dt1 > 1e-9, dt1, np.nan)
-        dang1 = np.abs(np.diff(col)) / dt1
-        best1 = int(np.nanargmax(dang1))
-        result[f'{axis}_vel1_max(deg/s)'] = float(dang1[best1])
-        result[f'{axis}_vel1_ts']         = float(ts[best1])
-
-        # windowed velocity  (frame i → i+window)
-        if n > window:
-            dtW   = ts[window:] - ts[:n - window]
-            dtW   = np.where(dtW > 1e-9, dtW, np.nan)
-            dangW = np.abs(col[window:] - col[:n - window]) / dtW
-            bestW = int(np.nanargmax(dangW))
-            result[f'{axis}_vel{window}_max(deg/s)'] = float(dangW[bestW])
-            result[f'{axis}_vel{window}_ts']         = float(ts[bestW])
+
+    # ── finite-diff rates (glitch-sensitive) ──────────────────────────────
+    dts = np.diff(ts)
+    dts = np.where(dts > 1e-9, dts, np.nan)
+    fd  = np.abs(np.diff(pose, axis=0)) / dts[:, np.newaxis]  # (N-1, 3) deg/s
+
+    for ci, ax in enumerate(axis_names):
+        result[f'{ax}_range(deg)'] = float(pose[:, ci].max() - pose[:, ci].min())
+        result[f'{ax}_fd_max(deg/s)'] = float(np.nanmax(fd[:, ci]))
+
+    # ── _dt rates (smooth, algorithm's own estimate) ──────────────────────
+    if dt_arr is not None and len(dt_arr) >= 2:
+        # align _dt timestamps to pose timestamps (inner join)
+        pose_ts_set = set(np.round(ts, 6))
+        dt_matched  = dt_arr[np.array([round(t, 6) in pose_ts_set
+                                       for t in dt_arr[:, 0]])]
+        if len(dt_matched) >= 2:
+            dt_ang = np.degrees(np.abs(dt_matched[:, 4:7]))  # deg/s, all positive
+            for ci, ax in enumerate(axis_names):
+                best = int(np.argmax(dt_ang[:, ci]))
+                result[f'{ax}_dt_max(deg/s)'] = float(dt_ang[best, ci])
+                result[f'{ax}_dt_ts']         = float(dt_matched[best, 0])
+        else:
+            for ax in axis_names:
+                result[f'{ax}_dt_max(deg/s)'] = float('nan')
+                result[f'{ax}_dt_ts']         = float('nan')
     else:
-            result[f'{axis}_vel{window}_max(deg/s)'] = float('nan')
-            result[f'{axis}_vel{window}_ts']         = float('nan')
+        for ax in axis_names:
+            result[f'{ax}_dt_max(deg/s)'] = float('nan')
+            result[f'{ax}_dt_ts']         = float('nan')
+
+    # ── glitch ratio: fd_max / dt_max ─────────────────────────────────────
+    # Values >> 1 indicate pose jumps not reflected in the smooth _dt fit.
+    # A ratio of ~5 matches a single inter-frame spike (vel5 ≈ vel1/5).
+    for ax in axis_names:
+        fd_max = result[f'{ax}_fd_max(deg/s)']
+        dt_max = result[f'{ax}_dt_max(deg/s)']
+        if dt_max and dt_max > 0.01:   # avoid divide-by-near-zero
+            result[f'{ax}_glitch'] = float(fd_max / dt_max)
+        else:
+            result[f'{ax}_glitch'] = float('nan')
 
     return result
 
@@ -166,7 +214,7 @@ def find_xml_files(root):
         if not os.path.isdir(scene_path):
             continue
 
-        # find version subdirs (e.g. v88, v80, …) — pick the latest by name
+        # find version subdirs (e.g. v88, v80) — pick latest by name
         try:
             subdirs = [d for d in os.listdir(scene_path)
                        if os.path.isdir(os.path.join(scene_path, d))
@@ -175,44 +223,40 @@ def find_xml_files(root):
             continue
         if not subdirs:
             continue
-        ver = sorted(subdirs)[-1]   # lexicographic; "v88" > "v80" > "v20" etc.
+        ver = sorted(subdirs)[-1]
 
         ver_path = os.path.join(scene_path, ver)
-        xml_name = f"{scene}-INTERFRAME.corr-xml"
-        xml_path = os.path.join(ver_path, xml_name)
+        xml_path = os.path.join(ver_path, f"{scene}-INTERFRAME.corr-xml")
         if os.path.isfile(xml_path):
             yield scene, xml_path
 
 
-# ──────────────────────────────── main ───────────────────────────────────── #
+# ─────────────────────────────── main ────────────────────────────────────── #
 
 def main():
     ap = argparse.ArgumentParser(description=__doc__,
                                  formatter_class=argparse.RawDescriptionHelpFormatter)
     ap.add_argument('--root',       required=True,
-                    help='Linked data root, e.g. .../linked_1763232117-1763234145-v88')
-    ap.add_argument('--window', type=int, default=5,
-                    help='Frame offset for windowed velocity (default 5)')
+                    help='Linked data root')
     ap.add_argument('--top',        type=int,   default=20,
-                    help='How many top sequences to print per axis (default 20)')
+                    help='Top sequences to print per axis (default 20)')
     ap.add_argument('--out',        default=None,
                     help='Output TSV path (optional)')
-    ap.add_argument('--sort',   default='R_vel1',
-                    choices=['R_range', 'R_vel1', 'R_velW',
-                             'T_range', 'T_vel1', 'T_velW',
-                             'A_range', 'A_vel1', 'A_velW'],
-                    help='Column to rank sequences by (default R_vel1)')
+    ap.add_argument('--sort',       default='R_dt',
+                    choices=['R_dt','T_dt','A_dt','R_range','T_range','A_range'],
+                    help='Column to rank by (default R_dt)')
+    ap.add_argument('--max_glitch', type=float, default=None,
+                    help='Hide sequences whose R glitch ratio exceeds this '
+                         '(e.g. 5.0 keeps only clean data; default: show all)')
     args = ap.parse_args()
 
-    W = args.window
     print(f"Scanning {args.root} …", flush=True)
-
     rows = []
     for scene_ts, xml_path in find_xml_files(args.root):
-        arr = parse_xml(xml_path)
-        if arr is None:
+        pose_arr, dt_arr = parse_xml(xml_path)
+        if pose_arr is None:
             continue
-        m = sequence_metrics(arr, W)
+        m = sequence_metrics(pose_arr, dt_arr)
         m['scene_ts'] = scene_ts
         m['xml_path'] = xml_path
         rows.append(m)
@@ -221,52 +265,57 @@ def main():
     if not rows:
         return
 
-    # ── sort map ──────────────────────────────────────────────────────────
-    sort_col_map = {
-        'R_range': f'r_range(deg)',
-        'R_vel1':  f'r_vel1_max(deg/s)',
-        'R_velW':  f'r_vel{W}_max(deg/s)',
-        'T_range': f't_range(deg)',
-        'T_vel1':  f't_vel1_max(deg/s)',
-        'T_velW':  f't_vel{W}_max(deg/s)',
-        'A_range': f'a_range(deg)',
-        'A_vel1':  f'a_vel1_max(deg/s)',
-        'A_velW':  f'a_vel{W}_max(deg/s)',
+    sort_map = {
+        'R_dt':    'r_dt_max(deg/s)',
+        'T_dt':    't_dt_max(deg/s)',
+        'A_dt':    'a_dt_max(deg/s)',
+        'R_range': 'r_range(deg)',
+        'T_range': 't_range(deg)',
+        'A_range': 'a_range(deg)',
     }
-    sort_key = sort_col_map[args.sort]
-    rows.sort(key=lambda r: r.get(sort_key, 0) or 0, reverse=True)
+    sort_key = sort_map[args.sort]
+
+    # apply glitch filter
+    display = rows
+    if args.max_glitch is not None:
+        display = [r for r in rows
+                   if (r.get('r_glitch') or 0) <= args.max_glitch]
+        print(f"After glitch filter (r_glitch ≤ {args.max_glitch}): "
+              f"{len(display)} sequences remain.")
+
+    display.sort(key=lambda r: r.get(sort_key) or 0, reverse=True)
 
     # ── console report ────────────────────────────────────────────────────
-    # For each axis print top-N by vel1 and by range
+    hdr = (f"  {'Scene':>22}  {'fr':>4}  {'range(°)':>8}  "
+           f"{'dt_max(°/s)':>11}  {'at timestamp':>20}  "
+           f"{'glitch':>7}  {'fd_max(°/s)':>11}")
+
     for axis, label in (('r', 'R (heading/course)'),
                         ('t', 'T (banking)'),
                         ('a', 'A (pitch)')):
-        print(f"\n{'═'*70}")
-        print(f"  {label}  —  top {args.top} by consecutive velocity")
-        print(f"{'═'*70}")
-        print(f"  {'Scene':>22}  {'frames':>6}  {'range(°)':>9}  "
-              f"{'vel1(°/s)':>10}  {'at timestamp':>20}  "
-              f"{'vel{W}(°/s)':>10}  {'at timestamp':>20}".replace('{W}', str(W)))
-        by_vel = sorted(rows, key=lambda r: r.get(f'{axis}_vel1_max(deg/s)', 0) or 0,
+        ranked = sorted(display,
+                        key=lambda r: r.get(f'{axis}_dt_max(deg/s)') or 0,
                         reverse=True)
-        for r in by_vel[:args.top]:
-            print(f"  {r['scene_ts']:>22}  {r['n_frames']:>6}  "
-                  f"{r[f'{axis}_range(deg)']:>9.2f}  "
-                  f"{r[f'{axis}_vel1_max(deg/s)']:>10.3f}  "
-                  f"{r[f'{axis}_vel1_ts']:>20.3f}  "
-                  f"{r[f'{axis}_vel{W}_max(deg/s)']:>10.3f}  "
-                  f"{r[f'{axis}_vel{W}_ts']:>20.3f}")
+        print(f"\n{'═'*90}")
+        print(f"  {label}  —  top {args.top} by _dt angular rate  "
+              f"(glitch = fd_max / dt_max; ~1 = clean)")
+        print(f"{'═'*90}")
+        print(hdr)
+        for r in ranked[:args.top]:
+            print(f"  {r['scene_ts']:>22}  {r['n_frames']:>4}  "
+                  f"{r[f'{axis}_range(deg)']:>8.2f}  "
+                  f"{r[f'{axis}_dt_max(deg/s)']:>11.3f}  "
+                  f"{r[f'{axis}_dt_ts']:>20.3f}  "
+                  f"{r[f'{axis}_glitch']:>7.1f}  "
+                  f"{r[f'{axis}_fd_max(deg/s)']:>11.1f}")
 
     # ── TSV output ────────────────────────────────────────────────────────
     if args.out:
         fieldnames = [
             'scene_ts', 'n_frames', 'duration_s',
-            'r_range(deg)', 'r_vel1_max(deg/s)', 'r_vel1_ts',
-            f'r_vel{W}_max(deg/s)', f'r_vel{W}_ts',
-            't_range(deg)', 't_vel1_max(deg/s)', 't_vel1_ts',
-            f't_vel{W}_max(deg/s)', f't_vel{W}_ts',
-            'a_range(deg)', 'a_vel1_max(deg/s)', 'a_vel1_ts',
-            f'a_vel{W}_max(deg/s)', f'a_vel{W}_ts',
+            'r_dt_max(deg/s)', 'r_dt_ts', 'r_range(deg)', 'r_glitch', 'r_fd_max(deg/s)',
+            't_dt_max(deg/s)', 't_dt_ts', 't_range(deg)', 't_glitch', 't_fd_max(deg/s)',
+            'a_dt_max(deg/s)', 'a_dt_ts', 'a_range(deg)', 'a_glitch', 'a_fd_max(deg/s)',
             'xml_path',
         ]
         sep = '\t' if args.out.endswith('.tsv') else ','
@@ -274,8 +323,10 @@ def main():
             w = csv.DictWriter(f, fieldnames=fieldnames, delimiter=sep,
                                extrasaction='ignore')
             w.writeheader()
-            w.writerows(rows)
-        print(f"\nSaved: {args.out}")
+            # write all rows (unfiltered) sorted by chosen key
+            for r in sorted(rows, key=lambda r: r.get(sort_key) or 0, reverse=True):
+                w.writerow(r)
+        print(f"\nSaved: {args.out}  ({len(rows)} rows, unfiltered)")
 
 
 if __name__ == '__main__':