CLAUDE: Add COLMAP sanity-check script

Verifies Phase 0 output before any OpenMVS tools run:
- camera intrinsics (focal length, FoV)
- camera-centre trajectory plots (XY/XZ/YZ)
- inter-frame step statistics with outlier detection
- consecutive rotation change (warns if > 10°)
- sample image display with NaN/saturation diagnostics

Fixed a stride-2 parse bug: empty 2D-point lines must be kept
in the line list before the range(0,N,2) loop, otherwise only
every other image record is parsed.
Co-authored-by: Claude <claude@elphel.com>

scripts/check_colmap.py

+#!/usr/bin/env python3
+# -----------------------------------------------------------------------------
+# check_colmap.py — Sanity-check a COLMAP sparse model from Elphel data
+#
+# Copyright (C) 2026 Elphel, Inc.
+# -----------------------------------------------------------------------------
+#
+# check_colmap.py is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+# -----------------------------------------------------------------------------
+"""
+Visualisation and sanity checks for the Phase 0 COLMAP output.
+
+What is verified
+----------------
+1. Camera-centre trajectory  (XY, XZ, YZ planes + 3-D scatter)
+   Expected: ~straight flight path, X-range ≈ 150 m over ~185 frames.
+
+2. Inter-frame step size
+   Expected: ~0.7 m between adjacent frames (60 Hz, ~44 m/s).
+   Outliers flag pose discontinuities.
+
+3. Rotation spread
+   Prints min/max/mean of pairwise angular distance between consecutive
+   frames.  Should be < 5° for this dataset.
+
+4. Sample images  (first / middle / last scene PNG)
+   Displays the 16-bit PNG slices.  Checks that they are not all-black
+   or all-white (a common sign of wrong normalisation).
+
+5. cameras.txt summary
+   Confirms focal length, image size, model type.
+
+Usage
+-----
+  python3 scripts/check_colmap.py --colmap <colmap_dir>
+
+  # e.g.
+  python3 scripts/check_colmap.py --colmap /tmp/colmap_test
+
+Dependencies: numpy, matplotlib, Pillow
+"""
+
+import argparse
+import sys
+from pathlib import Path
+
+import numpy as np
+import matplotlib
+matplotlib.use("Agg")   # headless by default; change to "TkAgg" for interactive
+import matplotlib.pyplot as plt
+from PIL import Image
+
+
+# ── COLMAP file readers ───────────────────────────────────────────────────────
+
+def read_cameras(cameras_txt):
+    cameras = {}
+    with open(cameras_txt) as f:
+        for line in f:
+            line = line.strip()
+            if not line or line.startswith("#"):
+                continue
+            parts = line.split()
+            cam_id = int(parts[0])
+            model  = parts[1]
+            w, h   = int(parts[2]), int(parts[3])
+            params = [float(x) for x in parts[4:]]
+            cameras[cam_id] = dict(model=model, width=w, height=h, params=params)
+    return cameras
+
+
+def _quat_to_R(qw, qx, qy, qz):
+    """Quaternion (w, x, y, z) → 3×3 rotation matrix (world-to-camera)."""
+    q = np.array([qw, qx, qy, qz])
+    q /= np.linalg.norm(q)
+    w, x, y, z = q
+    return np.array([
+        [1 - 2*(y*y + z*z),     2*(x*y - w*z),     2*(x*z + w*y)],
+        [    2*(x*y + w*z), 1 - 2*(x*x + z*z),     2*(y*z - w*x)],
+        [    2*(x*z - w*y),     2*(y*z + w*x), 1 - 2*(x*x + y*y)],
+    ])
+
+
+def read_images(images_txt):
+    """
+    Returns list of dicts, one per image, in file order:
+      id, qw, qx, qy, qz, tx, ty, tz, camera_id, name, R, centre
+    R      : 3×3 world-to-camera rotation matrix
+    centre : camera centre in world coords  C = −R^T @ t
+    """
+    entries = []
+    with open(images_txt) as f:
+        lines = [l.strip() for l in f if not l.startswith("#")]
+    # images.txt has two lines per image; even lines are data, odd lines are
+    # the (empty) 2D-point list.  Keep blank lines so the stride-2 pattern works.
+    for i in range(0, len(lines), 2):
+        if not lines[i]:
+            continue
+        parts = lines[i].split()
+        img_id   = int(parts[0])
+        qw, qx, qy, qz = float(parts[1]), float(parts[2]), float(parts[3]), float(parts[4])
+        tx, ty, tz      = float(parts[5]), float(parts[6]), float(parts[7])
+        cam_id   = int(parts[8])
+        name     = parts[9]
+        R = _quat_to_R(qw, qx, qy, qz)
+        t = np.array([tx, ty, tz])
+        centre = -R.T @ t
+        entries.append(dict(
+            id=img_id, qw=qw, qx=qx, qy=qy, qz=qz,
+            tx=tx, ty=ty, tz=tz, camera_id=cam_id, name=name,
+            R=R, centre=centre,
+        ))
+    return entries
+
+
+# ── Checks ────────────────────────────────────────────────────────────────────
+
+def check_cameras(cameras):
+    print("\n── cameras.txt ──────────────────────────────────────────────────")
+    for cam_id, c in cameras.items():
+        print(f"  Camera {cam_id}: {c['model']}  {c['width']}×{c['height']}")
+        params = c['params']
+        if c['model'] == 'PINHOLE':
+            fx, fy, cx, cy = params
+            print(f"    fx={fx:.2f}  fy={fy:.2f}  cx={cx:.2f}  cy={cy:.2f}")
+            hfov = 2 * np.degrees(np.arctan2(c['width'] / 2, fx))
+            vfov = 2 * np.degrees(np.arctan2(c['height'] / 2, fy))
+            print(f"    HFoV={hfov:.1f}°  VFoV={vfov:.1f}°")
+        else:
+            print(f"    params: {params}")
+
+
+def check_trajectory(entries, out_dir):
+    centres = np.array([e['centre'] for e in entries])   # (N, 3)
+
+    x, y, z = centres[:, 0], centres[:, 1], centres[:, 2]
+    print("\n── Camera centres (world frame) ─────────────────────────────────")
+    print(f"  N frames : {len(entries)}")
+    print(f"  X  {x.min():.2f} … {x.max():.2f} m   range {x.max()-x.min():.2f} m")
+    print(f"  Y  {y.min():.2f} … {y.max():.2f} m   range {y.max()-y.min():.2f} m")
+    print(f"  Z  {z.min():.2f} … {z.max():.2f} m   range {z.max()-z.min():.2f} m")
+
+    # Inter-frame step sizes
+    steps = np.linalg.norm(np.diff(centres, axis=0), axis=1)
+    print(f"\n── Inter-frame step (metres) ─────────────────────────────────────")
+    print(f"  mean={steps.mean():.3f}  median={np.median(steps):.3f}"
+          f"  min={steps.min():.3f}  max={steps.max():.3f}")
+    n_outliers = np.sum(steps > 3 * np.median(steps))
+    if n_outliers:
+        print(f"  WARNING: {n_outliers} steps > 3× median — check pose gaps")
+    else:
+        print("  No outlier steps detected.")
+
+    # Trajectory plots
+    fig, axes = plt.subplots(1, 3, figsize=(15, 4))
+    fig.suptitle("Camera trajectory (world frame)")
+
+    def _arrow(ax, xi, yi, stride=10):
+        for i in range(0, len(xi) - stride, stride):
+            ax.annotate("", xy=(xi[i+stride], yi[i+stride]),
+                        xytext=(xi[i], yi[i]),
+                        arrowprops=dict(arrowstyle="->", color="gray", lw=0.7))
+
+    axes[0].plot(x, y, "b.-", ms=2, lw=0.5); _arrow(axes[0], x, y)
+    axes[0].set_xlabel("X (m)"); axes[0].set_ylabel("Y (m)"); axes[0].set_title("XY (top view)")
+    axes[0].set_aspect("equal"); axes[0].grid(True, lw=0.3)
+
+    axes[1].plot(x, z, "b.-", ms=2, lw=0.5); _arrow(axes[1], x, z)
+    axes[1].set_xlabel("X (m)"); axes[1].set_ylabel("Z (m)"); axes[1].set_title("XZ (side view)")
+    axes[1].grid(True, lw=0.3)
+
+    axes[2].plot(y, z, "b.-", ms=2, lw=0.5); _arrow(axes[2], y, z)
+    axes[2].set_xlabel("Y (m)"); axes[2].set_ylabel("Z (m)"); axes[2].set_title("YZ (cross-track)")
+    axes[2].grid(True, lw=0.3)
+
+    # Mark first / reference / last
+    for ax, xi, yi in [(axes[0], x, y), (axes[1], x, z), (axes[2], y, z)]:
+        ax.plot(xi[0],  yi[0],  "go", ms=6, label="first")
+        ax.plot(xi[-1], yi[-1], "rs", ms=6, label="last")
+        ax.legend(fontsize=7)
+
+    plt.tight_layout()
+    out = out_dir / "trajectory.png"
+    plt.savefig(out, dpi=150)
+    plt.close()
+    print(f"\n  Trajectory plot → {out}")
+
+
+def check_rotations(entries):
+    print("\n── Consecutive rotation change (degrees) ────────────────────────")
+    angles = []
+    for i in range(1, len(entries)):
+        R1 = entries[i-1]['R']
+        R2 = entries[i  ]['R']
+        dR = R2 @ R1.T
+        # angle from rotation matrix: cos(θ) = (trace - 1) / 2
+        cos_angle = np.clip((np.trace(dR) - 1.0) / 2.0, -1.0, 1.0)
+        angles.append(np.degrees(np.arccos(cos_angle)))
+    angles = np.array(angles)
+    print(f"  mean={angles.mean():.3f}°  max={angles.max():.3f}°  "
+          f"99th-pct={np.percentile(angles, 99):.3f}°")
+    if angles.max() > 10:
+        print("  WARNING: large rotation steps detected — check Euler angle convention")
+
+
+def check_sample_images(entries, images_dir, out_dir, n_samples=3):
+    print(f"\n── Sample images ({n_samples} frames) ───────────────────────────────────")
+    n = len(entries)
+    indices = [0, n // 2, n - 1]
+    fig, axes = plt.subplots(1, n_samples, figsize=(5 * n_samples, 4))
+    for ax, idx in zip(axes, indices):
+        e = entries[idx]
+        img_path = images_dir / e['name']
+        if not img_path.exists():
+            ax.set_title(f"MISSING: {e['name']}")
+            ax.axis("off")
+            continue
+        img = np.array(Image.open(img_path), dtype=np.float32)
+        vmin, vmax = np.nanpercentile(img, 1), np.nanpercentile(img, 99)
+        ax.imshow(img, cmap="gray", vmin=vmin, vmax=vmax)
+        ax.set_title(f"Frame {idx+1}: {e['name']}\n"
+                     f"range [{img.min():.0f}, {img.max():.0f}]  "
+                     f"1-99%: [{vmin:.0f}, {vmax:.0f}]", fontsize=7)
+        ax.axis("off")
+        # Black / white saturation warning
+        black_frac = np.mean(img == 0)
+        white_frac = np.mean(img == 65535)
+        # 90% threshold: end-frames of terrain data normally have high NaN→0 fill
+        if black_frac > 0.9:
+            print(f"  WARNING frame {idx+1}: {black_frac*100:.0f}% pixels are 0 "
+                  "(image may be under-exposed or NaN-flooded)")
+        else:
+            print(f"  INFO frame {idx+1}: {black_frac*100:.0f}% pixels are 0 "
+                  "(NaN-filled regions)")
+        if white_frac > 0.1:
+            print(f"  WARNING frame {idx+1}: {white_frac*100:.0f}% pixels saturate "
+                  "at 65535")
+
+    plt.tight_layout()
+    out = out_dir / "sample_images.png"
+    plt.savefig(out, dpi=150)
+    plt.close()
+    print(f"  Sample image grid → {out}")
+
+
+# ── Main ─────────────────────────────────────────────────────────────────────
+
+def main():
+    ap = argparse.ArgumentParser(
+        description="Sanity-check COLMAP output from elphel_to_colmap.py")
+    ap.add_argument("--colmap", required=True,
+                    help="COLMAP output directory (contains images/ and sparse/0/)")
+    ap.add_argument("--out", default=None,
+                    help="Directory for output plots (default: <colmap>/checks/)")
+    args = ap.parse_args()
+
+    colmap_dir = Path(args.colmap)
+    sparse_dir = colmap_dir / "sparse" / "0"
+    images_dir = colmap_dir / "images"
+    out_dir    = Path(args.out) if args.out else colmap_dir / "checks"
+    out_dir.mkdir(parents=True, exist_ok=True)
+
+    for p in [sparse_dir / "cameras.txt", sparse_dir / "images.txt", images_dir]:
+        if not p.exists():
+            print(f"ERROR: expected path not found: {p}", file=sys.stderr)
+            sys.exit(1)
+
+    cameras = read_cameras(sparse_dir / "cameras.txt")
+    entries = read_images(sparse_dir / "images.txt")
+    print(f"Loaded {len(entries)} images from {sparse_dir / 'images.txt'}")
+
+    check_cameras(cameras)
+    check_trajectory(entries, out_dir)
+    check_rotations(entries)
+    check_sample_images(entries, images_dir, out_dir)
+
+    print(f"\nAll checks done.  Plots written to {out_dir}/")
+
+
+if __name__ == "__main__":
+    main()