Merge branch 'lane/b-recon'

backend/festival4d/frames.py

@@ -11,18 +11,78 @@ from __future__ import annotations
 import logging
 from pathlib import Path
 
+import numpy as np
+
 log = logging.getLogger("festival4d.frames")
 
 
 def sharpness(image) -> float:
-    """Variance of the Laplacian of a grayscale image (OpenCV) — higher is sharper."""
+    """Variance of the Laplacian of an image (OpenCV) — higher is sharper.
-    raise NotImplementedError("lane B (M2): implement sharpness")
+
+    Accepts a BGR or grayscale ``ndarray`` (as returned by ``cv2.VideoCapture.read``).
+    The variance of the Laplacian is the standard focus/blur measure: a sharp frame has
+    strong high-frequency edges (high variance), a blurred one is smooth (low variance).
+    """
+    import cv2
+
+    img = np.asarray(image)
+    if img.ndim == 3:
+        img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
+    return float(cv2.Laplacian(img, cv2.CV_64F).var())
 
 
 def sample_frames(video_path: Path, video_id: int, out_dir: Path,
                   target_fps: float = 2.0, window_s: float = 0.5) -> list[Path]:
     """Sample the sharpest frame per ``window_s`` window from ``video_path`` (spec M2).
 
-    Returns the written JPEG paths (named ``{video_id}_{frame_idx}.jpg``).
+    Decodes every frame, buckets them into non-overlapping ``window_s``-second windows, and
+    keeps the single sharpest frame (variance of Laplacian) in each window. This yields
+    ~``1 / window_s`` frames per second (≈ ``target_fps`` at the defaults), biased toward the
+    in-focus frames COLMAP needs. Writes JPEGs named ``{video_id}_{frame_idx}.jpg`` (the
+    original decoded ``frame_idx``, so ``t_video = frame_idx / fps``) into ``out_dir``.
+
+    Returns the written JPEG paths, ordered by frame index.
     """
-    raise NotImplementedError("lane B (M2): implement sample_frames")
+    import cv2
+
+    out_dir = Path(out_dir)
+    out_dir.mkdir(parents=True, exist_ok=True)
+
+    cap = cv2.VideoCapture(str(video_path))
+    if not cap.isOpened():
+        raise RuntimeError(f"could not open video for frame sampling: {video_path}")
+    try:
+        fps = float(cap.get(cv2.CAP_PROP_FPS))
+        if not np.isfinite(fps) or fps <= 0.0:
+            log.warning("frames: %s reported fps=%s; falling back to 30.0", video_path, fps)
+            fps = 30.0
+        bucket_s = float(window_s) if window_s and window_s > 0 else 0.5
+
+        # window_idx -> (sharpness, frame_idx, frame_bgr). Only the current best per window
+        # is retained, so memory stays ~O(number of windows), not O(frames).
+        best: dict[int, tuple[float, int, np.ndarray]] = {}
+        frame_idx = 0
+        while True:
+            ok, frame = cap.read()
+            if not ok:
+                break
+            t = frame_idx / fps
+            win = int(t / bucket_s)
+            s = sharpness(frame)
+            cur = best.get(win)
+            if cur is None or s > cur[0]:
+                best[win] = (s, frame_idx, frame)
+            frame_idx += 1
+    finally:
+        cap.release()
+
+    written: list[Path] = []
+    for win in sorted(best):
+        _, fidx, frame = best[win]
+        path = out_dir / f"{video_id}_{fidx}.jpg"
+        if not cv2.imwrite(str(path), frame):
+            raise RuntimeError(f"failed to write frame JPEG: {path}")
+        written.append(path)
+    log.info("frames: %s -> %d frames sampled into %s", Path(video_path).name,
+             len(written), out_dir)
+    return written

backend/festival4d/geometry.py

@@ -121,7 +121,8 @@ def colmap_to_threejs(q: ArrayLike, t: ArrayLike) -> tuple[NDArray[np.float64],
 
 
 # ---------------------------------------------------------------------------
-# Lane B stubs (spec M2 + M8). Signatures FROZEN; bodies raise NotImplementedError.
+# Lane B (spec M2 + M8): pose interpolation, pixel ray casting, two-view triangulation,
+# and the single-view nearest-point fallback. Signatures FROZEN; bodies implemented below.
 # ---------------------------------------------------------------------------
 def slerp_pose(
     q0: ArrayLike,
@@ -149,7 +150,39 @@ def slerp_pose(
     -------
     (q, t) : the interpolated quaternion ``[w, x, y, z]`` and translation ``[x, y, z]``.
     """
-    raise NotImplementedError("lane B (M2): implement slerp_pose")
+    a = float(alpha)
+    q0 = np.asarray(q0, dtype=np.float64).reshape(4)
+    q1 = np.asarray(q1, dtype=np.float64).reshape(4)
+    n0 = np.linalg.norm(q0)
+    n1 = np.linalg.norm(q1)
+    if n0 < 1e-12 or n1 < 1e-12:
+        raise ValueError("slerp endpoint quaternion has near-zero norm")
+    q0 = q0 / n0
+    q1 = q1 / n1
+
+    # Double cover: pick the sign of q1 that lies on the same hemisphere as q0, so slerp
+    # takes the shorter arc (a rotation and its negation are the same orientation).
+    dot = float(np.dot(q0, q1))
+    if dot < 0.0:
+        q1 = -q1
+        dot = -dot
+    dot = min(1.0, max(-1.0, dot))
+
+    if dot > 0.9995:
+        # Endpoints almost coincide: nlerp is numerically safe and visually identical.
+        q_interp = q0 + a * (q1 - q0)
+    else:
+        theta_0 = np.arccos(dot)
+        sin_0 = np.sin(theta_0)
+        s0 = np.sin((1.0 - a) * theta_0) / sin_0
+        s1 = np.sin(a * theta_0) / sin_0
+        q_interp = s0 * q0 + s1 * q1
+    q_interp = q_interp / np.linalg.norm(q_interp)
+
+    t0 = np.asarray(t0, dtype=np.float64).reshape(3)
+    t1 = np.asarray(t1, dtype=np.float64).reshape(3)
+    t_interp = (1.0 - a) * t0 + a * t1
+    return q_interp, t_interp
 
 
 def ray_from_pixel(
@@ -174,7 +207,20 @@ def ray_from_pixel(
     direction : ndarray, shape (3,)
         Unit ray direction in world coords, pointing into the scene.
     """
-    raise NotImplementedError("lane B (M8): implement ray_from_pixel")
+    R = quat_to_mat(q)                          # world -> cam
+    t_vec = np.asarray(t, dtype=np.float64).reshape(3)
+    origin = -R.T @ t_vec                        # camera center in world coords
+
+    # Pinhole back-projection. A world point X projects with x_cam = R X + t and
+    # px = fx * x_cam.x / x_cam.z + cx, py = fy * x_cam.y / x_cam.z + cy
+    # (COLMAP camera axes: +x right, +y down, +z forward). So the camera-space direction
+    # through pixel (px, py) is [(px-cx)/fx, (py-cy)/fy, 1], pointing forward into the scene.
+    d_cam = np.array([(px - cx) / fx, (py - cy) / fy, 1.0], dtype=np.float64)
+    d_world = R.T @ d_cam                         # rotate direction cam -> world
+    norm = np.linalg.norm(d_world)
+    if norm < 1e-12:
+        raise ValueError("degenerate ray direction")
+    return origin, d_world / norm
 
 
 def triangulate_rays(
@@ -194,7 +240,37 @@ def triangulate_rays(
         reject the triangulation when the rays are near-parallel or ``gap`` exceeds the
         spec threshold (0.5 scene units).
     """
-    raise NotImplementedError("lane B (M8): implement triangulate_rays")
+    oa = np.asarray(origin_a, dtype=np.float64).reshape(3)
+    ob = np.asarray(origin_b, dtype=np.float64).reshape(3)
+    da = np.asarray(dir_a, dtype=np.float64).reshape(3)
+    db = np.asarray(dir_b, dtype=np.float64).reshape(3)
+    na, nb = np.linalg.norm(da), np.linalg.norm(db)
+    if na < 1e-12 or nb < 1e-12:
+        raise ValueError("triangulate_rays: zero-length direction")
+    da = da / na
+    db = db / nb
+
+    # Shortest segment between two lines P(s)=oa+s*da, Q(u)=ob+u*db. Minimize |P-Q|^2.
+    # With unit directions: b = da.db, denom = 1 - b^2 (0 when parallel).
+    w0 = oa - ob
+    b = float(np.dot(da, db))
+    d = float(np.dot(da, w0))
+    e = float(np.dot(db, w0))
+    denom = 1.0 - b * b
+    if denom < 1e-9:
+        # Near-parallel: no unique closest pair. Anchor on oa, take the closest point on
+        # line b to it; gap is the line-to-line perpendicular distance. Callers reject
+        # near-parallel rays up front, so this branch just stays numerically safe.
+        s = 0.0
+        u = e
+    else:
+        s = (b * e - d) / denom
+        u = (e - b * d) / denom
+    pa = oa + s * da
+    pb = ob + u * db
+    point = 0.5 * (pa + pb)
+    gap = float(np.linalg.norm(pa - pb))
+    return point, gap
 
 
 def nearest_point_on_ray(
@@ -214,4 +290,26 @@ def nearest_point_on_ray(
     point : ndarray shape (3,) or None
         The selected point-cloud point, or ``None`` if none lie within ``radius``.
     """
-    raise NotImplementedError("lane B (M8): implement nearest_point_on_ray")
+    o = np.asarray(origin, dtype=np.float64).reshape(3)
+    d = np.asarray(direction, dtype=np.float64).reshape(3)
+    nd = np.linalg.norm(d)
+    if nd < 1e-12:
+        raise ValueError("nearest_point_on_ray: zero-length direction")
+    d = d / nd
+
+    pts = np.asarray(points, dtype=np.float64).reshape(-1, 3)
+    if len(pts) == 0:
+        return None
+
+    v = pts - o                          # origin -> each point
+    proj = v @ d                         # signed distance along the ray
+    perp = v - np.outer(proj, d)         # component perpendicular to the ray
+    perp_dist = np.linalg.norm(perp, axis=1)
+
+    # In front of the origin and inside the cylinder of the given radius.
+    mask = (proj > 0.0) & (perp_dist <= radius)
+    if not np.any(mask):
+        return None
+    idx_in = np.where(mask)[0]
+    best = idx_in[np.argmin(perp_dist[idx_in])]
+    return pts[best].copy()

backend/festival4d/sfm.py

@@ -1,33 +1,61 @@
 """COLMAP orchestration, model parsing, scene normalization, pose export (spec M2, lane B).
 
-STUB — lane B fills the bodies; the public signatures below are frozen (``cli.py`` calls
+Pipeline (spec M2):
-``run_reconstruct``).
-
-Pipeline (spec M2, for lane B's reference):
 - Sample frames (``frames.sample_frames``), then drive the COLMAP CLI via subprocess:
-  ``feature_extractor`` (OPENCV model, single camera per folder) -> ``sequential_matcher``
+  ``feature_extractor`` (OPENCV model, single camera per folder) -> ``exhaustive_matcher``
-  within each video + ``exhaustive_matcher`` across -> ``mapper`` -> ``model_converter`` to
+  (small frame counts; exhaustive covers within- and cross-video pairs) -> ``mapper`` ->
-  a TXT model. Parse ``images.txt`` / ``cameras.txt`` / ``points3D.txt`` by hand (no
+  ``image_undistorter`` + ``model_converter`` to a PINHOLE TXT model. Parse ``images.txt`` /
-  pycolmap).
+  ``cameras.txt`` / ``points3D.txt`` by hand (no pycolmap).
 - **Normalize** the scene: centroid -> origin, camera bounding sphere radius -> 10, world-up
-  -> average camera up. Apply the same similarity transform to poses and points.
+  -> average camera up. The same similarity transform is applied to poses and points.
 - **Interpolate** poses for unregistered frames (slerp rotation, lerp translation, mark
-  ``registered=False``); never extrapolate past the first/last registered frame.
+  ``registered=False``); never extrapolate past the first/last registered frame of a video.
-- Export ``config.POINTS_PLY`` (binary little-endian, matching ``synthetic.write_ply``'s
+- Export ``config.POINTS_PLY`` (binary little-endian, matching ``synthetic.write_ply``) and
-  format) and fill ``camera_poses`` via ``db.set_poses``.
+  fill ``camera_poses`` via ``db.set_poses``.
-- **Failure handling** (spec M2): COLMAP optional at runtime (``shutil.which("colmap")``);
+- **Failure handling** (spec M2): COLMAP optional at runtime (``shutil.which("colmap")``); if
-  if it registers < 60% of frames or < 2 videos into one model, print a clear diagnostic and
+  it registers < 60% of sampled frames or < 2 videos into one model, print a clear diagnostic
-  leave existing poses untouched — never corrupt the DB. The app degrades to "synced videos,
+  and leave existing poses untouched — never corrupt the DB. The app degrades to "synced
-  no 3D".
+  videos, no 3D".
 """
 
 from __future__ import annotations
 
 import logging
+import functools
+import struct
+import subprocess
 from pathlib import Path
 
+import numpy as np
+
+from festival4d import config, db, frames
+from festival4d.geometry import mat_to_quat, quat_to_mat, slerp_pose
+
 log = logging.getLogger("festival4d.sfm")
 
+# Minimum share of sampled frames COLMAP must register, and minimum videos in one model,
+# for the reconstruction to be trusted (spec M2 failure handling).
+MIN_REGISTERED_FRACTION = 0.60
+MIN_VIDEOS_IN_MODEL = 2
+NORMALIZED_SPHERE_RADIUS = 10.0
+
+# How COLMAP camera models lay out PARAMS -> (fx, fy, cx, cy). Distortion terms are ignored:
+# the stored intrinsics are the pinhole part (exact after image_undistorter's PINHOLE output,
+# approximate if a distorted model is parsed directly). Covers the models COLMAP emits here.
+_INTRINSICS_FROM_PARAMS = {
+    "SIMPLE_PINHOLE": lambda p: (p[0], p[0], p[1], p[2]),
+    "PINHOLE": lambda p: (p[0], p[1], p[2], p[3]),
+    "SIMPLE_RADIAL": lambda p: (p[0], p[0], p[1], p[2]),
+    "RADIAL": lambda p: (p[0], p[0], p[1], p[2]),
+    "SIMPLE_RADIAL_FISHEYE": lambda p: (p[0], p[0], p[1], p[2]),
+    "RADIAL_FISHEYE": lambda p: (p[0], p[0], p[1], p[2]),
+    "OPENCV": lambda p: (p[0], p[1], p[2], p[3]),
+    "OPENCV_FISHEYE": lambda p: (p[0], p[1], p[2], p[3]),
+    "FULL_OPENCV": lambda p: (p[0], p[1], p[2], p[3]),
+    "FOV": lambda p: (p[0], p[1], p[2], p[3]),
+    "THIN_PRISM_FISHEYE": lambda p: (p[0], p[1], p[2], p[3]),
+}
+
 
 def colmap_available() -> bool:
     """Whether the COLMAP binary is on PATH (spec M2 optional-at-runtime rule)."""
@@ -36,29 +64,416 @@ def colmap_available() -> bool:
     return shutil.which("colmap") is not None
 
 
+# ---------------------------------------------------------------------------
+# COLMAP TXT-model parsers (hand-written; no pycolmap dependency)
+# ---------------------------------------------------------------------------
+def _iter_data_lines(path: Path):
+    """Yield non-empty, non-comment lines from a COLMAP TXT file."""
+    with open(path, "r", encoding="utf-8") as f:
+        for line in f:
+            line = line.strip()
+            if line and not line.startswith("#"):
+                yield line
+
+
 def parse_images_txt(path: Path) -> list[dict]:
-    """Parse a COLMAP ``images.txt`` into per-image pose records (spec M2, lane B)."""
+    """Parse a COLMAP ``images.txt`` into per-image pose records (spec M2, lane B).
-    raise NotImplementedError("lane B (M2): implement parse_images_txt")
+
+    Format: two lines per image. The first is
+    ``IMAGE_ID QW QX QY QZ TX TY TZ CAMERA_ID NAME``; the second is the 2D keypoint list
+    (ignored here). Returns one dict per image with the world->camera quaternion
+    ``[w,x,y,z]``, translation, camera id, and image name.
+    """
+    records: list[dict] = []
+    take_pose_line = True  # image data alternates: pose line, then points2D line
+    for line in _iter_data_lines(path):
+        if take_pose_line:
+            parts = line.split()
+            if len(parts) < 10:
+                raise ValueError(f"malformed images.txt pose line: {line!r}")
+            records.append({
+                "image_id": int(parts[0]),
+                "qw": float(parts[1]), "qx": float(parts[2]),
+                "qy": float(parts[3]), "qz": float(parts[4]),
+                "tx": float(parts[5]), "ty": float(parts[6]), "tz": float(parts[7]),
+                "camera_id": int(parts[8]),
+                "name": " ".join(parts[9:]),  # NAME may itself contain spaces
+            })
+        take_pose_line = not take_pose_line
+    return records
 
 
 def parse_cameras_txt(path: Path) -> dict[int, dict]:
-    """Parse a COLMAP ``cameras.txt`` into ``camera_id -> intrinsics`` (spec M2, lane B)."""
+    """Parse a COLMAP ``cameras.txt`` into ``camera_id -> intrinsics`` (spec M2, lane B).
-    raise NotImplementedError("lane B (M2): implement parse_cameras_txt")
+
+    Format per line: ``CAMERA_ID MODEL WIDTH HEIGHT PARAMS[]``. Returns, per camera id, a
+    dict with ``model``, ``width``, ``height``, pinhole ``fx, fy, cx, cy`` and the raw
+    ``params`` list.
+    """
+    cameras: dict[int, dict] = {}
+    for line in _iter_data_lines(path):
+        parts = line.split()
+        if len(parts) < 4:
+            raise ValueError(f"malformed cameras.txt line: {line!r}")
+        cam_id = int(parts[0])
+        model = parts[1]
+        width, height = int(parts[2]), int(parts[3])
+        params = [float(x) for x in parts[4:]]
+        if model not in _INTRINSICS_FROM_PARAMS:
+            raise ValueError(f"unsupported COLMAP camera model {model!r}")
+        try:
+            fx, fy, cx, cy = _INTRINSICS_FROM_PARAMS[model](params)
+        except IndexError:
+            raise ValueError(
+                f"cameras.txt: model {model} has too few PARAMS: {line!r}"
+            ) from None
+        cameras[cam_id] = {
+            "model": model, "width": width, "height": height,
+            "fx": float(fx), "fy": float(fy), "cx": float(cx), "cy": float(cy),
+            "params": params,
+        }
+    return cameras
 
 
-def parse_points3d_txt(path: Path):
+def parse_points3d_txt(path: Path) -> tuple[np.ndarray, np.ndarray]:
-    """Parse a COLMAP ``points3D.txt`` into ``(points Nx3, colors Nx3)`` (spec M2, lane B)."""
+    """Parse a COLMAP ``points3D.txt`` into ``(points Nx3 float32, colors Nx3 uint8)``.
-    raise NotImplementedError("lane B (M2): implement parse_points3d_txt")
+
+    Format per line: ``POINT3D_ID X Y Z R G B ERROR TRACK[]`` (track ignored).
+    """
+    xyz: list[tuple[float, float, float]] = []
+    rgb: list[tuple[int, int, int]] = []
+    for line in _iter_data_lines(path):
+        parts = line.split()
+        if len(parts) < 7:
+            raise ValueError(f"malformed points3D.txt line: {line!r}")
+        xyz.append((float(parts[1]), float(parts[2]), float(parts[3])))
+        rgb.append((int(parts[4]), int(parts[5]), int(parts[6])))
+    if not xyz:
+        return (np.zeros((0, 3), dtype=np.float32), np.zeros((0, 3), dtype=np.uint8))
+    return (np.asarray(xyz, dtype=np.float32), np.asarray(rgb, dtype=np.uint8))
+
+
+# ---------------------------------------------------------------------------
+# Scene normalization
+# ---------------------------------------------------------------------------
+def _rotation_aligning(a: np.ndarray, b: np.ndarray) -> np.ndarray:
+    """Shortest-arc rotation matrix that maps unit vector ``a`` onto unit vector ``b``."""
+    a = a / np.linalg.norm(a)
+    b = b / np.linalg.norm(b)
+    v = np.cross(a, b)
+    c = float(np.dot(a, b))
+    s = float(np.linalg.norm(v))
+    if s < 1e-12:
+        if c > 0:                       # already aligned
+            return np.eye(3)
+        # antiparallel: rotate 180 deg about any axis perpendicular to a
+        perp = np.array([1.0, 0.0, 0.0])
+        if abs(a[0]) > 0.9:
+            perp = np.array([0.0, 1.0, 0.0])
+        axis = np.cross(a, perp)
+        axis /= np.linalg.norm(axis)
+        return 2.0 * np.outer(axis, axis) - np.eye(3)
+    vx = np.array([[0.0, -v[2], v[1]], [v[2], 0.0, -v[0]], [-v[1], v[0], 0.0]])
+    return np.eye(3) + vx + vx @ vx * ((1.0 - c) / (s * s))
+
+
+def _camera_center(pose: dict) -> np.ndarray:
+    """World-space camera center ``C = -R^T t`` for a world->camera pose dict."""
+    R = quat_to_mat([pose["qw"], pose["qx"], pose["qy"], pose["qz"]])
+    t = np.array([pose["tx"], pose["ty"], pose["tz"]], dtype=np.float64)
+    return -R.T @ t
+
+
+def _camera_up(pose: dict) -> np.ndarray:
+    """World-space camera up vector. COLMAP camera +y is *down*, so up = R^T @ (0,-1,0)."""
+    R = quat_to_mat([pose["qw"], pose["qx"], pose["qy"], pose["qz"]])
+    return R.T @ np.array([0.0, -1.0, 0.0])
 
 
 def normalize_scene(points, poses):
-    """Similarity transform: centroid->origin, camera sphere radius->10, up->+Y (spec M2)."""
+    """Similarity transform: centroid->origin, camera sphere radius->10, up->+Y (spec M2).
-    raise NotImplementedError("lane B (M2): implement normalize_scene")
+
+    ``points`` is an ``(N,3)`` array; ``poses`` a list of world->camera pose dicts (keys
+    ``qw,qx,qy,qz,tx,ty,tz`` plus any extras, which are preserved). The single similarity
+    transform ``X' = s * Rn @ (X - c)`` (translate to the point centroid, rotate the average
+    camera-up onto +Y, scale so the camera bounding sphere has radius 10) is applied to both
+    the points and the poses.
+
+    Returns ``(points_new float32 (N,3), poses_new list[dict])``.
+    """
+    pts = np.asarray(points, dtype=np.float64).reshape(-1, 3)
+    poses = list(poses)
+
+    # 1. Translate to the point-cloud centroid.
+    centroid = pts.mean(axis=0) if len(pts) else np.zeros(3)
+
+    # 2. Rotate so the average camera-up aligns with +Y (phones held roughly upright).
+    if poses:
+        ups = np.array([_camera_up(p) for p in poses])
+        up_avg = ups.mean(axis=0)
+        if np.linalg.norm(up_avg) > 1e-9:
+            Rn = _rotation_aligning(up_avg, np.array([0.0, 1.0, 0.0]))
+        else:
+            Rn = np.eye(3)
+    else:
+        Rn = np.eye(3)
+
+    # 3. Scale so the farthest camera center from the centroid sits at radius 10.
+    if poses:
+        centers = np.array([_camera_center(p) for p in poses])
+        max_dist = float(np.max(np.linalg.norm(centers - centroid, axis=1)))
+        scale = NORMALIZED_SPHERE_RADIUS / max_dist if max_dist > 1e-9 else 1.0
+    else:
+        scale = 1.0
+
+    # Apply X' = s * Rn @ (X - c) to points.
+    points_new = (scale * (pts - centroid) @ Rn.T).astype(np.float32)
+
+    # Apply the matching transform to each world->camera pose:
+    #   R' = R @ Rn^T ,  t' = s * (t + R @ c)   (derived so projections are unchanged).
+    poses_new: list[dict] = []
+    for p in poses:
+        R = quat_to_mat([p["qw"], p["qx"], p["qy"], p["qz"]])
+        t = np.array([p["tx"], p["ty"], p["tz"]], dtype=np.float64)
+        R_new = R @ Rn.T
+        t_new = scale * (t + R @ centroid)
+        q_new = mat_to_quat(R_new)
+        out = dict(p)
+        out["qw"], out["qx"], out["qy"], out["qz"] = (
+            float(q_new[0]), float(q_new[1]), float(q_new[2]), float(q_new[3]))
+        out["tx"], out["ty"], out["tz"] = float(t_new[0]), float(t_new[1]), float(t_new[2])
+        poses_new.append(out)
+
+    return points_new, poses_new
 
 
-def run_colmap(frames_dir: Path, workspace: Path) -> Path:
+# ---------------------------------------------------------------------------
-    """Drive the COLMAP CLI and return the exported TXT-model directory (spec M2, lane B)."""
+# Pose interpolation
-    raise NotImplementedError("lane B (M2): implement run_colmap")
+# ---------------------------------------------------------------------------
+def interpolate_poses(registered: list[dict], sampled: list[dict]) -> list[dict]:
+    """Fill unregistered sampled frames of one video between registered ones (spec M2).
+
+    ``registered`` are solved poses (with ``registered=True``); ``sampled`` are the frames
+    that were fed to COLMAP as ``{frame_idx, t_video_s, ...}``. For each sampled frame that
+    lies strictly between the first and last registered frame but wasn't solved, slerp the
+    rotation and lerp the translation of the two bracketing registered poses (intrinsics
+    inherited from the earlier neighbor), marking it ``registered=False``. Frames outside the
+    registered span are dropped (no extrapolation). Returns all poses sorted by ``t_video_s``.
+    """
+    reg = sorted(registered, key=lambda p: p["t_video_s"])
+    if not reg:
+        return []
+    reg_frames = {int(p["frame_idx"]) for p in reg}
+    t_first, t_last = reg[0]["t_video_s"], reg[-1]["t_video_s"]
+
+    out: list[dict] = list(reg)
+    for sf in sampled:
+        fidx = int(sf["frame_idx"])
+        tv = float(sf["t_video_s"])
+        if fidx in reg_frames:
+            continue
+        if tv <= t_first or tv >= t_last:
+            continue  # no extrapolation beyond the registered span
+        # bracketing registered neighbors
+        lo = max((p for p in reg if p["t_video_s"] <= tv), key=lambda p: p["t_video_s"])
+        hi = min((p for p in reg if p["t_video_s"] >= tv), key=lambda p: p["t_video_s"])
+        span = hi["t_video_s"] - lo["t_video_s"]
+        alpha = 0.0 if span <= 1e-12 else (tv - lo["t_video_s"]) / span
+        q, t = slerp_pose(
+            [lo["qw"], lo["qx"], lo["qy"], lo["qz"]], [lo["tx"], lo["ty"], lo["tz"]],
+            [hi["qw"], hi["qx"], hi["qy"], hi["qz"]], [hi["tx"], hi["ty"], hi["tz"]],
+            alpha,
+        )
+        out.append({
+            "frame_idx": fidx, "t_video_s": tv,
+            "qw": float(q[0]), "qx": float(q[1]), "qy": float(q[2]), "qz": float(q[3]),
+            "tx": float(t[0]), "ty": float(t[1]), "tz": float(t[2]),
+            "fx": lo["fx"], "fy": lo["fy"], "cx": lo["cx"], "cy": lo["cy"],
+            "registered": False,
+        })
+    out.sort(key=lambda p: p["t_video_s"])
+    return out
+
+
+# ---------------------------------------------------------------------------
+# COLMAP CLI orchestration
+# ---------------------------------------------------------------------------
+@functools.lru_cache(maxsize=None)
+def _colmap_help(command: str) -> str:
+    """Cached ``colmap <command> --help`` text, for build-aware option detection.
+
+    COLMAP prints its option list to stderr, so both streams are captured.
+    """
+    try:
+        proc = subprocess.run(
+            ["colmap", command, "--help"], capture_output=True, text=True
+        )
+        return (proc.stdout or "") + "\n" + (proc.stderr or "")
+    except Exception:
+        return ""
+
+
+def _cpu_flag(command: str) -> list[str]:
+    """Force CPU SIFT for ``command`` if this COLMAP build exposes a ``use_gpu`` option.
+
+    The option prefix differs by version: COLMAP 3.x uses ``SiftExtraction`` /
+    ``SiftMatching``, 4.x uses ``FeatureExtraction`` / ``FeatureMatching``. We read the option
+    name from ``--help`` so the pipeline runs headless on either (and on CPU-only builds that
+    omit the option entirely, we pass nothing). Returns e.g. ``["--FeatureExtraction.use_gpu",
+    "0"]`` or ``[]``.
+    """
+    for line in _colmap_help(command).splitlines():
+        line = line.strip()
+        if line.startswith("--") and ".use_gpu" in line:
+            return [line.split()[0], "0"]
+    return []
+
+
+def _run_colmap_step(args: list[str]) -> bool:
+    """Run one ``colmap <args>`` step; return True on success, False (logged) on failure."""
+    cmd = ["colmap", *args]
+    log.info("colmap: %s", " ".join(args[:2]))
+    try:
+        proc = subprocess.run(cmd, capture_output=True, text=True)
+    except FileNotFoundError:
+        log.error("colmap binary vanished from PATH")
+        return False
+    if proc.returncode != 0:
+        log.error("colmap %s failed (rc=%d): %s", args[0], proc.returncode,
+                  (proc.stderr or proc.stdout or "").strip()[-500:])
+        return False
+    return True
+
+
+def _registered_image_count(model_dir: Path) -> int | None:
+    """Number of registered images in a COLMAP model dir, or None if it has no image list.
+
+    The mapper writes binary models by default; ``images.bin`` begins with a little-endian
+    ``uint64`` giving the registered-image count, which is the exact metric we want (file
+    size would instead track total keypoint observations and mis-rank multi-component runs).
+    Falls back to counting the TXT model's two-lines-per-image list.
+    """
+    images_txt = model_dir / "images.txt"
+    images_bin = model_dir / "images.bin"
+    if images_txt.exists():
+        return sum(1 for _ in _iter_data_lines(images_txt)) // 2
+    if images_bin.exists():
+        with open(images_bin, "rb") as f:
+            header = f.read(8)
+        return struct.unpack("<Q", header)[0] if len(header) == 8 else 0
+    return None
+
+
+def _largest_model_dir(sparse_dir: Path) -> Path | None:
+    """Return the reconstruction subdir (``0``, ``1``, ...) with the most registered images."""
+    candidates = [d for d in sparse_dir.iterdir() if d.is_dir()] if sparse_dir.exists() else []
+    best, best_n = None, -1
+    for d in candidates:
+        n = _registered_image_count(d)
+        if n is None:
+            continue
+        if n > best_n:
+            best, best_n = d, n
+    return best
+
+
+def run_colmap(frames_dir: Path, workspace: Path) -> Path | None:
+    """Drive the COLMAP CLI and return the exported TXT-model directory (spec M2, lane B).
+
+    ``frames_dir`` holds one subfolder of JPEGs per video (single camera per folder). Runs
+    feature extraction (CPU SIFT, OPENCV model), exhaustive matching, incremental mapping,
+    then undistortion + TXT conversion. Returns the directory containing ``images.txt`` /
+    ``cameras.txt`` / ``points3D.txt``, or ``None`` if COLMAP produced no usable model
+    (caller degrades gracefully).
+    """
+    workspace = Path(workspace)
+    # Start each run from a clean workspace: a stale database.db makes feature_extractor
+    # fail ("images already exist"), and a leftover sparse/ model would poison the pick.
+    if workspace.exists():
+        import shutil
+        shutil.rmtree(workspace)
+    workspace.mkdir(parents=True, exist_ok=True)
+    database = workspace / "database.db"
+    sparse = workspace / "sparse"
+    sparse.mkdir(exist_ok=True)
+
+    ok = _run_colmap_step([
+        "feature_extractor",
+        "--database_path", str(database),
+        "--image_path", str(frames_dir),
+        "--ImageReader.single_camera_per_folder", "1",
+        "--ImageReader.camera_model", "OPENCV",
+        *_cpu_flag("feature_extractor"),
+    ])
+    if not ok:
+        return None
+
+    # Exhaustive matching covers within- and cross-video pairs; frame counts are small.
+    if not _run_colmap_step([
+        "exhaustive_matcher",
+        "--database_path", str(database),
+        *_cpu_flag("exhaustive_matcher"),
+    ]):
+        return None
+
+    if not _run_colmap_step([
+        "mapper",
+        "--database_path", str(database),
+        "--image_path", str(frames_dir),
+        "--output_path", str(sparse),
+    ]):
+        return None
+
+    model = _largest_model_dir(sparse)
+    if model is None:
+        log.error("colmap mapper produced no reconstruction")
+        return None
+
+    # Undistort -> PINHOLE model, then convert to TXT. Fall back to converting the raw
+    # (possibly distorted) sparse model directly if undistortion fails.
+    txt_dir = workspace / "model_txt"
+    txt_dir.mkdir(exist_ok=True)
+    dense = workspace / "dense"
+    undistorted = _run_colmap_step([
+        "image_undistorter",
+        "--image_path", str(frames_dir),
+        "--input_path", str(model),
+        "--output_path", str(dense),
+        "--output_type", "COLMAP",
+    ])
+    convert_input = (dense / "sparse") if undistorted and (dense / "sparse").exists() else model
+    if not _run_colmap_step([
+        "model_converter",
+        "--input_path", str(convert_input),
+        "--output_path", str(txt_dir),
+        "--output_type", "TXT",
+    ]):
+        return None
+    if not (txt_dir / "images.txt").exists():
+        log.error("colmap model_converter produced no images.txt")
+        return None
+    return txt_dir
+
+
+# ---------------------------------------------------------------------------
+# Full reconstruction entrypoint
+# ---------------------------------------------------------------------------
+def _video_id_and_frame(name: str) -> tuple[int, int] | None:
+    """Recover ``(video_id, frame_idx)`` from a ``{video_id}_{frame_idx}.jpg`` image name."""
+    stem = Path(name).stem
+    if "_" not in stem:
+        return None
+    vid_s, frame_s = stem.rsplit("_", 1)
+    try:
+        return int(vid_s), int(frame_s)
+    except ValueError:
+        return None
+
+
+def _diagnostic(msg: str) -> None:
+    """Emit a clear operator-facing diagnostic (spec M2 graceful degradation)."""
+    log.warning("reconstruct: %s", msg)
+    print(f"[reconstruct] {msg}")
 
 
 def run_reconstruct() -> dict:
@@ -66,6 +481,128 @@ def run_reconstruct() -> dict:
 
     Entrypoint for ``python -m festival4d reconstruct``. Degrades gracefully when COLMAP is
     absent or the reconstruction is too weak (spec M2 failure handling): existing poses are
-    left untouched. Returns a summary dict.
+    left untouched and the DB is never corrupted. Returns a summary dict with a ``status``.
     """
-    raise NotImplementedError("lane B (M2): implement run_reconstruct")
+    db.init_engine()
+    db.init_db()
+    videos = db.get_videos()
+    if not videos:
+        _diagnostic("no videos in the project — run `synthetic` or `ingest` first.")
+        return {"status": "no_videos"}
+
+    if not colmap_available():
+        _diagnostic(
+            "COLMAP not found on PATH — skipping 3D reconstruction. Install it "
+            "(`brew install colmap` on macOS) to enable pose/point-cloud solving. "
+            "Existing (synthetic/previous) poses are left untouched; the app still runs as "
+            "synced videos without 3D."
+        )
+        return {"status": "skipped_no_colmap", "videos": len(videos)}
+
+    # 1. Sample sharp frames per video into one subfolder each (single camera per folder).
+    frames_root = config.FRAMES_DIR
+    per_video_sampled: dict[int, list[dict]] = {}
+    for v in videos:
+        video_path = config.RAW_DIR / v.filename
+        if not video_path.exists():
+            _diagnostic(f"video file missing, skipping: {video_path}")
+            continue
+        out_dir = frames_root / str(v.id)
+        paths = frames.sample_frames(video_path, v.id, out_dir)
+        fps = float(v.fps) if v.fps and v.fps > 0 else 30.0
+        sampled = []
+        for p in paths:
+            parsed = _video_id_and_frame(p.name)
+            if parsed is None:
+                continue
+            _, fidx = parsed
+            sampled.append({"frame_idx": fidx, "t_video_s": fidx / fps})
+        per_video_sampled[v.id] = sorted(sampled, key=lambda s: s["frame_idx"])
+
+    total_sampled = sum(len(s) for s in per_video_sampled.values())
+    if total_sampled == 0:
+        _diagnostic("no frames could be sampled from any video — nothing to reconstruct.")
+        return {"status": "no_frames"}
+
+    # 2. Run COLMAP.
+    model_dir = run_colmap(frames_root, config.COLMAP_DIR)
+    if model_dir is None:
+        _diagnostic(
+            "COLMAP produced no usable reconstruction. Likely causes: too-dark footage, "
+            "motion blur, or insufficient view overlap between cameras. Existing poses left "
+            "untouched — the app degrades to synced videos, no 3D."
+        )
+        return {"status": "failed_no_model", "sampled": total_sampled}
+
+    # 3. Parse the model.
+    images = parse_images_txt(model_dir / "images.txt")
+    cameras = parse_cameras_txt(model_dir / "cameras.txt")
+    points, colors = parse_points3d_txt(model_dir / "points3D.txt")
+
+    # 4. Attach each registered image to its video + intrinsics.
+    registered_by_video: dict[int, list[dict]] = {}
+    for img in images:
+        parsed = _video_id_and_frame(img["name"])
+        if parsed is None:
+            continue
+        video_id, frame_idx = parsed
+        cam = cameras.get(img["camera_id"])
+        if cam is None:
+            continue
+        fps = next((float(v.fps) for v in videos if v.id == video_id), 30.0) or 30.0
+        registered_by_video.setdefault(video_id, []).append({
+            "frame_idx": frame_idx, "t_video_s": frame_idx / fps,
+            "qw": img["qw"], "qx": img["qx"], "qy": img["qy"], "qz": img["qz"],
+            "tx": img["tx"], "ty": img["ty"], "tz": img["tz"],
+            "fx": cam["fx"], "fy": cam["fy"], "cx": cam["cx"], "cy": cam["cy"],
+            "registered": True, "video_id": video_id,
+        })
+
+    registered_count = sum(len(v) for v in registered_by_video.values())
+    videos_in_model = len(registered_by_video)
+    frac = registered_count / total_sampled if total_sampled else 0.0
+    if frac < MIN_REGISTERED_FRACTION or videos_in_model < MIN_VIDEOS_IN_MODEL:
+        _diagnostic(
+            f"weak reconstruction: {registered_count}/{total_sampled} frames "
+            f"({frac:.0%}) registered across {videos_in_model} video(s); need "
+            f">= {MIN_REGISTERED_FRACTION:.0%} of frames and >= {MIN_VIDEOS_IN_MODEL} "
+            "videos. Likely causes: too-dark footage, motion blur, or too little view "
+            "overlap. Existing poses left untouched — degrading to synced videos, no 3D."
+        )
+        return {
+            "status": "failed_weak", "sampled": total_sampled,
+            "registered": registered_count, "videos_in_model": videos_in_model,
+        }
+
+    # 5. Normalize the whole scene (points + all registered poses) with one transform.
+    all_registered = [p for v in registered_by_video.values() for p in v]
+    norm_points, norm_registered = normalize_scene(points, all_registered)
+    norm_by_video: dict[int, list[dict]] = {}
+    for p in norm_registered:
+        norm_by_video.setdefault(p["video_id"], []).append(p)
+
+    # 6. Interpolate unregistered sampled frames, then write poses per video (atomic replace).
+    interpolated_total = 0
+    for video_id, reg in norm_by_video.items():
+        sampled = per_video_sampled.get(video_id, [])
+        final = interpolate_poses(reg, sampled)
+        interpolated_total += sum(1 for p in final if not p.get("registered", True))
+        db.set_poses(video_id, final)
+
+    # 7. Export the normalized point cloud in the frozen PLY format.
+    from festival4d.synthetic import write_ply
+    write_ply(config.POINTS_PLY, norm_points, colors)
+
+    summary = {
+        "status": "ok",
+        "videos_in_model": videos_in_model,
+        "registered": registered_count,
+        "interpolated": interpolated_total,
+        "points": int(len(norm_points)),
+        "points_ply": str(config.POINTS_PLY),
+    }
+    log.info("reconstruct: done — %s", summary)
+    print(f"[reconstruct] reconstructed {videos_in_model} videos, "
+          f"{registered_count} registered + {interpolated_total} interpolated poses, "
+          f"{len(norm_points)} points -> {config.POINTS_PLY}")
+    return summary

backend/tests/test_geometry.py

@@ -108,3 +108,226 @@ def test_quat_to_mat_matches_scipy_oracle():
         np.testing.assert_allclose(
             quat_to_mat(q_wxyz), Rotation.from_quat(q_xyzw).as_matrix(), atol=1e-9
         )
+
+
+# ===========================================================================
+# Lane B (spec M2 + M8) tests — below the frozen block. These exercise the
+# stubs foundation left: slerp_pose, ray_from_pixel, triangulate_rays,
+# nearest_point_on_ray. They do NOT touch the frozen colmap_to_threejs cases above.
+# ===========================================================================
+from scipy.spatial.transform import Slerp  # noqa: E402
+
+from festival4d.geometry import (  # noqa: E402
+    colmap_to_threejs,
+    nearest_point_on_ray,
+    ray_from_pixel,
+    slerp_pose,
+    triangulate_rays,
+)
+
+
+def _wxyz_to_xyzw(q):
+    q = np.asarray(q, dtype=float)
+    return np.array([q[1], q[2], q[3], q[0]])
+
+
+def _rand_unit_quat(rng):
+    q = rng.normal(size=4)
+    return q / np.linalg.norm(q)
+
+
+# --- slerp_pose ------------------------------------------------------------
+def test_slerp_pose_endpoints():
+    q0 = _rand_unit_quat(np.random.default_rng(1))
+    q1 = _rand_unit_quat(np.random.default_rng(2))
+    t0 = np.array([1.0, -2.0, 3.0])
+    t1 = np.array([-4.0, 5.0, 6.0])
+
+    q_a, t_a = slerp_pose(q0, t0, q1, t1, 0.0)
+    q_b, t_b = slerp_pose(q0, t0, q1, t1, 1.0)
+    # endpoints recover the endpoint *rotations* (quaternion up to sign) and translations
+    np.testing.assert_allclose(quat_to_mat(q_a), quat_to_mat(q0), atol=1e-12)
+    np.testing.assert_allclose(quat_to_mat(q_b), quat_to_mat(q1), atol=1e-12)
+    np.testing.assert_allclose(t_a, t0, atol=1e-12)
+    np.testing.assert_allclose(t_b, t1, atol=1e-12)
+
+
+def test_slerp_pose_matches_scipy_oracle():
+    rng = np.random.default_rng(20260716)
+    for _ in range(30):
+        q0 = _rand_unit_quat(rng)
+        q1 = _rand_unit_quat(rng)
+        oracle = Slerp([0.0, 1.0], Rotation.from_quat(
+            [_wxyz_to_xyzw(q0), _wxyz_to_xyzw(q1)]))
+        for alpha in (0.1, 0.25, 0.5, 0.73, 0.9):
+            q, _ = slerp_pose(q0, [0, 0, 0], q1, [0, 0, 0], alpha)
+            np.testing.assert_allclose(
+                quat_to_mat(q), oracle(alpha).as_matrix(), atol=1e-9)
+
+
+def test_slerp_pose_translation_is_linear():
+    q = np.array([1.0, 0.0, 0.0, 0.0])
+    t0 = np.array([0.0, 0.0, 0.0])
+    t1 = np.array([10.0, -4.0, 2.0])
+    for alpha in (0.0, 0.3, 0.5, 1.0):
+        _, t = slerp_pose(q, t0, q, t1, alpha)
+        np.testing.assert_allclose(t, (1 - alpha) * t0 + alpha * t1, atol=1e-12)
+
+
+def test_slerp_pose_double_cover_takes_short_arc():
+    """q1 and -q1 are the same rotation; slerp must yield the same (short-arc) result."""
+    rng = np.random.default_rng(99)
+    q0 = _rand_unit_quat(rng)
+    q1 = _rand_unit_quat(rng)
+    for alpha in (0.2, 0.5, 0.8):
+        qa, _ = slerp_pose(q0, [0, 0, 0], q1, [0, 0, 0], alpha)
+        qb, _ = slerp_pose(q0, [0, 0, 0], -q1, [0, 0, 0], alpha)
+        np.testing.assert_allclose(quat_to_mat(qa), quat_to_mat(qb), atol=1e-12)
+
+
+def test_slerp_pose_midpoint_is_half_angle():
+    """A 180-degree-ish pair: the midpoint rotation angle is half the endpoint angle."""
+    q0 = np.array([1.0, 0.0, 0.0, 0.0])                 # identity
+    ang = np.radians(100.0)
+    q1 = np.array([np.cos(ang / 2), 0.0, np.sin(ang / 2), 0.0])  # yaw 100 deg about Y
+    qm, _ = slerp_pose(q0, [0, 0, 0], q1, [0, 0, 0], 0.5)
+    # rotation angle of qm relative to identity should be ~50 deg
+    angle = 2.0 * np.arccos(min(1.0, abs(qm[0])))
+    assert abs(np.degrees(angle) - 50.0) < 1e-6
+
+
+# --- ray_from_pixel --------------------------------------------------------
+def test_ray_from_pixel_origin_is_camera_center():
+    """The ray origin equals the Three.js camera center from the frozen contract."""
+    rng = np.random.default_rng(5)
+    for _ in range(20):
+        q = _rand_unit_quat(rng)
+        t = rng.uniform(-5, 5, size=3)
+        position, _ = colmap_to_threejs(q, t)
+        origin, _ = ray_from_pixel(q, t, 600, 600, 320, 180, 320, 180)
+        np.testing.assert_allclose(origin, position, atol=1e-9)
+
+
+def test_ray_from_pixel_center_is_forward_axis():
+    """The center pixel unprojects along the camera forward (+z) axis, in world coords."""
+    q = np.array([1.0, 0.0, 0.0, 0.0])   # identity world->cam
+    t = np.array([0.0, 0.0, -10.0])      # camera center at (0,0,10)
+    origin, direction = ray_from_pixel(q, t, 500, 500, 320, 180, 320, 180)
+    np.testing.assert_allclose(origin, [0.0, 0.0, 10.0], atol=1e-9)
+    np.testing.assert_allclose(direction, [0.0, 0.0, 1.0], atol=1e-9)
+
+
+def _project_colmap(q, t, fx, fy, cx, cy, X):
+    """Forward COLMAP pinhole projection of a world point X -> pixel (px, py)."""
+    R = quat_to_mat(q)
+    xc = R @ np.asarray(X, float) + np.asarray(t, float)
+    return fx * xc[0] / xc[2] + cx, fy * xc[1] / xc[2] + cy, xc[2]
+
+
+def test_ray_from_pixel_is_projection_inverse():
+    """A world point in front of the camera lies exactly on the ray through its pixel."""
+    rng = np.random.default_rng(7)
+    for _ in range(50):
+        q = _rand_unit_quat(rng)
+        t = rng.uniform(-3, 3, size=3)
+        fx = fy = rng.uniform(400, 800)
+        cx, cy = 320.0, 180.0
+        C = -quat_to_mat(q).T @ t                    # camera center
+        forward = quat_to_mat(q).T @ np.array([0, 0, 1.0])
+        X = C + rng.uniform(2, 8) * forward + rng.uniform(-1, 1, size=3)  # in front
+        px, py, zc = _project_colmap(q, t, fx, fy, cx, cy, X)
+        if zc <= 0.1:
+            continue
+        origin, direction = ray_from_pixel(q, t, fx, fy, cx, cy, px, py)
+        to_X = X - origin
+        # X - origin must be parallel to direction and in front (positive projection)
+        cross = np.cross(to_X, direction)
+        assert np.linalg.norm(cross) < 1e-6 * (1 + np.linalg.norm(to_X))
+        assert np.dot(to_X, direction) > 0
+
+
+# --- triangulate_rays ------------------------------------------------------
+def test_triangulate_rays_intersecting():
+    P = np.array([1.0, 2.0, 3.0])
+    oa = np.array([0.0, 0.0, 0.0])
+    ob = np.array([4.0, 0.0, 0.0])
+    point, gap = triangulate_rays(oa, P - oa, ob, P - ob)
+    np.testing.assert_allclose(point, P, atol=1e-9)
+    assert gap < 1e-9
+
+
+def test_triangulate_rays_skew_known_geometry():
+    """Ray A along +x at z=0; ray B along +y at z=1. Closest points (0,0,0),(0,0,1)."""
+    oa = np.array([0.0, 0.0, 0.0]); da = np.array([1.0, 0.0, 0.0])
+    ob = np.array([0.0, 0.0, 1.0]); db = np.array([0.0, 1.0, 0.0])
+    point, gap = triangulate_rays(oa, da, ob, db)
+    np.testing.assert_allclose(point, [0.0, 0.0, 0.5], atol=1e-9)
+    assert abs(gap - 1.0) < 1e-9
+
+
+def test_triangulate_rays_parallel_is_safe():
+    """Parallel rays must not divide-by-zero; gap ~ their separation."""
+    oa = np.array([0.0, 0.0, 0.0]); da = np.array([1.0, 0.0, 0.0])
+    ob = np.array([0.0, 2.0, 0.0]); db = np.array([1.0, 0.0, 0.0])
+    point, gap = triangulate_rays(oa, da, ob, db)
+    assert np.all(np.isfinite(point))
+    assert abs(gap - 2.0) < 1e-6
+
+
+def test_triangulate_rays_from_two_cameras_recovers_point():
+    """Two cameras looking at a world point; triangulating their pixel rays recovers it."""
+    rng = np.random.default_rng(123)
+    P = np.array([0.5, 1.0, -0.5])
+    fx = fy = 600.0; cx, cy = 320.0, 180.0
+    rays = []
+    for center in ([6.0, 2.0, 6.0], [-6.0, 2.5, 6.0]):
+        C = np.array(center, float)
+        z = P - C; z /= np.linalg.norm(z)
+        up = np.array([0.0, 1.0, 0.0])
+        up = up - np.dot(up, z) * z; up /= np.linalg.norm(up)
+        x = np.cross(z, up); y = -up
+        R_c2w = np.column_stack([x, y, z]); R = R_c2w.T
+        t = -R @ C
+        q = mat_to_quat(R)
+        px, py, _ = _project_colmap(q, t, fx, fy, cx, cy, P)
+        rays.append(ray_from_pixel(q, t, fx, fy, cx, cy, px, py))
+    point, gap = triangulate_rays(rays[0][0], rays[0][1], rays[1][0], rays[1][1])
+    assert np.linalg.norm(point - P) < 0.2   # spec M8 acceptance tolerance
+    assert gap < 1e-6
+
+
+# --- nearest_point_on_ray --------------------------------------------------
+def test_nearest_point_on_ray_hits_within_radius():
+    o = np.array([0.0, 0.0, 0.0]); d = np.array([0.0, 0.0, 1.0])
+    points = np.array([
+        [0.1, 0.0, 5.0],    # perp 0.1, in front -> candidate
+        [0.05, 0.0, 3.0],   # perp 0.05, in front -> closest
+        [2.0, 0.0, 5.0],    # perp 2.0 -> outside radius
+    ])
+    got = nearest_point_on_ray(o, d, points, radius=0.3)
+    np.testing.assert_allclose(got, [0.05, 0.0, 3.0], atol=1e-12)
+
+
+def test_nearest_point_on_ray_excludes_outside_and_behind():
+    o = np.array([0.0, 0.0, 0.0]); d = np.array([0.0, 0.0, 1.0])
+    behind = np.array([[0.05, 0.0, -3.0]])     # in cylinder but behind origin
+    outside = np.array([[1.0, 0.0, 3.0]])      # in front but outside radius
+    assert nearest_point_on_ray(o, d, behind, radius=0.3) is None
+    assert nearest_point_on_ray(o, d, outside, radius=0.3) is None
+    assert nearest_point_on_ray(o, d, np.zeros((0, 3)), radius=0.3) is None
+
+
+def test_nearest_point_on_ray_synthetic_stage_corner():
+    """M8 single-view fallback: a ray toward a stage corner selects that corner point."""
+    from festival4d import synthetic
+
+    corner = np.array(synthetic.STAGE_CORNERS["Stage FL"], float)
+    points, _ = synthetic.generate_point_cloud()
+    # a camera looking straight at the corner; ray through the image center hits it
+    q, t = synthetic.look_at_colmap((5.0, 3.0, 9.0), corner)
+    intr = synthetic.intrinsics(640, 360)
+    origin, direction = ray_from_pixel(
+        q, t, intr["fx"], intr["fy"], intr["cx"], intr["cy"], intr["cx"], intr["cy"])
+    got = nearest_point_on_ray(origin, direction, points.astype(float), radius=0.3)
+    assert got is not None
+    np.testing.assert_allclose(got, corner, atol=1e-4)

backend/tests/test_sfm.py

@@ -0,0 +1,494 @@
+"""Tests for lane B: frames.py (sampling) + sfm.py (parsers, normalization,
+interpolation, and the graceful-degradation reconstruct pipeline).
+
+Everything here is verifiable without COLMAP: parsers run on hand-written TXT snippets, the
+export path runs on a hand-built (synthetic-pose) COLMAP model injected in place of a real
+COLMAP run, and the failure paths are exercised by monkeypatching. Spec M2 acceptance:
+"on the synthetic fixture (which skips COLMAP and injects known poses) the whole export path
+runs; if COLMAP is installed, the pipeline runs end-to-end without crashing."
+"""
+
+from __future__ import annotations
+
+import struct
+from pathlib import Path
+
+import numpy as np
+import pytest
+
+from festival4d import config, db, frames, sfm, synthetic
+from festival4d.geometry import quat_to_mat
+
+
+# ===========================================================================
+# COLMAP TXT parsers (hand-written snippets)
+# ===========================================================================
+IMAGES_TXT = """\
+# Image list with two lines of data per image:
+#   IMAGE_ID, QW, QX, QY, QZ, TX, TY, TZ, CAMERA_ID, NAME
+#   POINTS2D[] as (X, Y, POINT3D_ID)
+# Number of images: 2, mean observations per image: 2
+1 0.9998 0.0 0.02 0.0 -1.5 0.3 8.0 1 1/1_0.jpg
+480.0 270.0 5 500.1 260.2 -1
+2 0.7071 0.0 0.7071 0.0 2.0 0.1 7.5 2 2/2_30.jpg
+100.0 120.0 -1 300.0 200.0 5
+"""
+
+CAMERAS_TXT = """\
+# Camera list with one line of data per camera:
+#   CAMERA_ID, MODEL, WIDTH, HEIGHT, PARAMS[]
+# Number of cameras: 2
+1 PINHOLE 640 360 500.0 500.0 320.0 180.0
+2 OPENCV 640 360 510.0 511.0 321.0 181.0 0.01 -0.02 0.0 0.0
+"""
+
+POINTS3D_TXT = """\
+# 3D point list with one line of data per point:
+#   POINT3D_ID, X, Y, Z, R, G, B, ERROR, TRACK[] as (IMAGE_ID, POINT2D_IDX)
+# Number of points: 2, mean track length: 2
+5 1.0 2.0 3.0 200 100 50 0.5 1 0 2 1
+9 -1.0 0.0 4.0 10 20 30 0.8 1 1
+"""
+
+
+def test_parse_images_txt(tmp_path):
+    p = tmp_path / "images.txt"
+    p.write_text(IMAGES_TXT)
+    recs = sfm.parse_images_txt(p)
+    assert len(recs) == 2
+    assert recs[0]["image_id"] == 1
+    assert recs[0]["camera_id"] == 1
+    assert recs[0]["name"] == "1/1_0.jpg"
+    np.testing.assert_allclose(
+        [recs[0]["qw"], recs[0]["qx"], recs[0]["qy"], recs[0]["qz"]],
+        [0.9998, 0.0, 0.02, 0.0])
+    np.testing.assert_allclose(
+        [recs[1]["tx"], recs[1]["ty"], recs[1]["tz"]], [2.0, 0.1, 7.5])
+
+
+def test_parse_cameras_txt(tmp_path):
+    p = tmp_path / "cameras.txt"
+    p.write_text(CAMERAS_TXT)
+    cams = sfm.parse_cameras_txt(p)
+    assert set(cams) == {1, 2}
+    assert cams[1]["model"] == "PINHOLE"
+    assert (cams[1]["fx"], cams[1]["fy"], cams[1]["cx"], cams[1]["cy"]) == \
+        (500.0, 500.0, 320.0, 180.0)
+    # OPENCV: first 4 params are fx, fy, cx, cy; distortion ignored
+    assert (cams[2]["fx"], cams[2]["fy"], cams[2]["cx"], cams[2]["cy"]) == \
+        (510.0, 511.0, 321.0, 181.0)
+
+
+def test_parse_points3d_txt(tmp_path):
+    p = tmp_path / "points3D.txt"
+    p.write_text(POINTS3D_TXT)
+    points, colors = sfm.parse_points3d_txt(p)
+    assert points.shape == (2, 3) and colors.shape == (2, 3)
+    np.testing.assert_allclose(points[0], [1.0, 2.0, 3.0])
+    assert tuple(colors[0]) == (200, 100, 50)
+    assert points.dtype == np.float32 and colors.dtype == np.uint8
+
+
+def test_parse_points3d_empty(tmp_path):
+    p = tmp_path / "points3D.txt"
+    p.write_text("# header only\n")
+    points, colors = sfm.parse_points3d_txt(p)
+    assert points.shape == (0, 3) and colors.shape == (0, 3)
+
+
+# ===========================================================================
+# Scene normalization
+# ===========================================================================
+def _all_synth_poses():
+    poses = []
+    for i in range(3):
+        poses.extend(synthetic.camera_track(i, 20.0, 30, 640, 360))
+    return poses
+
+
+def test_normalize_scene_invariants():
+    points, _ = synthetic.generate_point_cloud()
+    poses = _all_synth_poses()
+    npoints, nposes = sfm.normalize_scene(points.astype(float), poses)
+
+    # centroid at origin
+    assert np.linalg.norm(npoints.mean(axis=0)) < 1e-3
+    # camera bounding sphere (from origin = point centroid) has radius 10
+    centers = np.array([sfm._camera_center(p) for p in nposes])
+    assert abs(float(np.max(np.linalg.norm(centers, axis=1))) - 10.0) < 1e-6
+    # average camera-up aligns with +Y
+    ups = np.array([sfm._camera_up(p) for p in nposes])
+    up_avg = ups.mean(axis=0)
+    up_avg /= np.linalg.norm(up_avg)
+    np.testing.assert_allclose(up_avg, [0.0, 1.0, 0.0], atol=1e-6)
+
+
+def test_normalize_scene_preserves_projection():
+    """A world point projects to the same pixel before and after normalization
+    (a similarity transform of the world must not change any image)."""
+    points, _ = synthetic.generate_point_cloud()
+    poses = _all_synth_poses()
+    intr = synthetic.intrinsics(640, 360)
+    for p in poses:
+        p.update(intr)
+    X = np.array([0.7, 0.9, 0.3])  # arbitrary world point
+
+    def project(pose, Xw):
+        R = quat_to_mat([pose["qw"], pose["qx"], pose["qy"], pose["qz"]])
+        t = np.array([pose["tx"], pose["ty"], pose["tz"]])
+        xc = R @ Xw + t
+        return np.array([pose["fx"] * xc[0] / xc[2] + pose["cx"],
+                         pose["fy"] * xc[1] / xc[2] + pose["cy"]])
+
+    npoints, nposes = sfm.normalize_scene(points.astype(float), poses)
+    # transform the same world point by X' = s*Rn@(X-c): recover s,Rn,c from the point map.
+    # Easier: the invariant is projection equality, so map X through the same transform used
+    # on the cloud by fitting it — instead, just check pixel equality using each pose pair.
+    # Reconstruct transform from three non-collinear cloud points is overkill; use the fact
+    # that projection is invariant, so compare original point's pixel to the transformed
+    # point's pixel where the transform is inferred from the point cloud centroid+scale.
+    c = points.astype(float).mean(axis=0)
+    centers = np.array([sfm._camera_center(p) for p in poses])
+    scale = 10.0 / float(np.max(np.linalg.norm(centers - c, axis=1)))
+    ups = np.array([sfm._camera_up(p) for p in poses])
+    up_avg = ups.mean(axis=0)
+    Rn = sfm._rotation_aligning(up_avg, np.array([0.0, 1.0, 0.0]))
+    Xn = scale * Rn @ (X - c)
+
+    for p_old, p_new in zip(poses, nposes):
+        np.testing.assert_allclose(project(p_old, X), project(p_new, Xn), atol=1e-6)
+
+
+# ===========================================================================
+# Pose interpolation
+# ===========================================================================
+def _pose(frame_idx, t, q, t3, registered=True):
+    return {"frame_idx": frame_idx, "t_video_s": t,
+            "qw": q[0], "qx": q[1], "qy": q[2], "qz": q[3],
+            "tx": t3[0], "ty": t3[1], "tz": t3[2],
+            "fx": 500.0, "fy": 500.0, "cx": 320.0, "cy": 180.0,
+            "registered": registered, "video_id": 1}
+
+
+def test_interpolate_poses_fills_gaps_no_extrapolation():
+    q0 = [1.0, 0.0, 0.0, 0.0]
+    ang = np.radians(60.0)
+    q1 = [np.cos(ang / 2), 0.0, np.sin(ang / 2), 0.0]  # 60 deg yaw about Y
+    registered = [_pose(0, 0.0, q0, [0, 0, 0]), _pose(60, 2.0, q1, [6, 0, 0])]
+    sampled = [{"frame_idx": f, "t_video_s": f / 30.0} for f in (0, 30, 60, 90)]
+
+    out = sfm.interpolate_poses(registered, sampled)
+    by_frame = {p["frame_idx"]: p for p in out}
+    assert set(by_frame) == {0, 30, 60}          # frame 90 dropped (beyond last registered)
+    assert by_frame[30]["registered"] is False    # interpolated
+    assert by_frame[0]["registered"] is True and by_frame[60]["registered"] is True
+    # midpoint translation is the lerp; rotation is ~30 deg
+    np.testing.assert_allclose(
+        [by_frame[30]["tx"], by_frame[30]["ty"], by_frame[30]["tz"]], [3, 0, 0], atol=1e-9)
+    mid_angle = 2.0 * np.degrees(np.arccos(min(1.0, abs(by_frame[30]["qw"]))))
+    assert abs(mid_angle - 30.0) < 1e-6
+
+
+def test_interpolate_poses_empty_registered():
+    assert sfm.interpolate_poses([], [{"frame_idx": 0, "t_video_s": 0.0}]) == []
+
+
+# ===========================================================================
+# frames.py — sharpness + windowed sampling
+# ===========================================================================
+def test_sharpness_sharp_beats_blurred():
+    import cv2
+
+    rng = np.random.default_rng(0)
+    noise = (rng.integers(0, 256, size=(120, 120), dtype=np.uint8))   # high-frequency
+    flat = np.full((120, 120), 128, dtype=np.uint8)                    # no edges
+    blurred = cv2.GaussianBlur(noise, (0, 0), sigmaX=5)
+    assert frames.sharpness(noise) > frames.sharpness(blurred) > frames.sharpness(flat)
+    # accepts BGR too
+    bgr = cv2.cvtColor(noise, cv2.COLOR_GRAY2BGR)
+    assert frames.sharpness(bgr) > 0
+
+
+def test_sample_frames_picks_sharpest_per_window(tmp_path):
+    import cv2
+
+    w = h = 128
+    fps = 30.0
+    n_frames = 30                       # 1.0 s -> two 0.5 s windows
+    vid_path = tmp_path / "clip.mp4"
+    writer = cv2.VideoWriter(str(vid_path), cv2.VideoWriter_fourcc(*"mp4v"), fps, (w, h))
+    if not writer.isOpened():
+        pytest.skip("no usable VideoWriter codec in this environment")
+    checker = np.indices((h, w)).sum(axis=0) % 2   # fine checkerboard -> sharp
+    sharp = (checker * 255).astype(np.uint8)
+    flat = np.full((h, w), 128, dtype=np.uint8)
+    sharp_frames = {7, 22}
+    for i in range(n_frames):
+        g = sharp if i in sharp_frames else flat
+        writer.write(cv2.cvtColor(g, cv2.COLOR_GRAY2BGR))
+    writer.release()
+
+    out = frames.sample_frames(vid_path, video_id=5, out_dir=tmp_path / "frames")
+    got_frames = sorted(int(Path(p).stem.split("_")[1]) for p in out)
+
+    # ground-truth: decode the written stream and take the per-window argmax ourselves
+    cap = cv2.VideoCapture(str(vid_path))
+    sh, idx = [], 0
+    while True:
+        ok, fr = cap.read()
+        if not ok:
+            break
+        sh.append(frames.sharpness(fr)); idx += 1
+    cap.release()
+    bucket = {}
+    for i, s in enumerate(sh):
+        win = int((i / fps) / 0.5)
+        if win not in bucket or s > bucket[win][0]:
+            bucket[win] = (s, i)
+    expected = sorted(v[1] for v in bucket.values())
+
+    assert got_frames == expected
+    assert all(Path(p).name.startswith("5_") for p in out)
+
+
+# ===========================================================================
+# run_reconstruct — graceful degradation + full export path (no real COLMAP)
+# ===========================================================================
+def _snapshot_poses():
+    return {
+        v.id: [(p.frame_idx, round(p.qw, 9), round(p.tx, 9), p.registered)
+               for p in db.get_poses(v.id)]
+        for v in db.get_videos()
+    }
+
+
+def _write_colmap_model(model_dir: Path, records: list[dict], cameras: dict,
+                        points, colors) -> None:
+    """Write a COLMAP TXT model (images/cameras/points3D) from pose records."""
+    model_dir.mkdir(parents=True, exist_ok=True)
+    with open(model_dir / "images.txt", "w") as f:
+        f.write("# image list\n")
+        for i, r in enumerate(records, start=1):
+            f.write(f"{i} {r['qw']} {r['qx']} {r['qy']} {r['qz']} "
+                    f"{r['tx']} {r['ty']} {r['tz']} {r['camera_id']} {r['name']}\n")
+            f.write("100.0 100.0 -1\n")  # dummy points2D line
+    with open(model_dir / "cameras.txt", "w") as f:
+        f.write("# camera list\n")
+        for cid, cam in cameras.items():
+            f.write(f"{cid} PINHOLE {cam['width']} {cam['height']} "
+                    f"{cam['fx']} {cam['fy']} {cam['cx']} {cam['cy']}\n")
+    with open(model_dir / "points3D.txt", "w") as f:
+        f.write("# point list\n")
+        for j, (p, c) in enumerate(zip(points, colors), start=1):
+            f.write(f"{j} {p[0]} {p[1]} {p[2]} {int(c[0])} {int(c[1])} {int(c[2])} 0.5 1 0\n")
+
+
+def test_reconstruct_no_videos():
+    db.init_engine()
+    db.reset_db()
+    result = sfm.run_reconstruct()
+    assert result["status"] == "no_videos"
+
+
+def test_reconstruct_skips_without_colmap(monkeypatch):
+    synthetic.build(run_ffmpeg=False)
+    monkeypatch.setattr(sfm, "colmap_available", lambda: False)
+    before = _snapshot_poses()
+    result = sfm.run_reconstruct()
+    assert result["status"] == "skipped_no_colmap"
+    assert _snapshot_poses() == before          # poses untouched
+
+
+def _fake_sampling(monkeypatch, frames_per_video):
+    """Make sampling hermetic: placeholder raw files + fake sample_frames output."""
+    for v in db.get_videos():
+        (config.RAW_DIR / v.filename).parent.mkdir(parents=True, exist_ok=True)
+        (config.RAW_DIR / v.filename).write_bytes(b"stub")
+
+    def fake_sample(video_path, video_id, out_dir, **kw):
+        return [Path(f"{video_id}_{f}.jpg") for f in frames_per_video]
+
+    monkeypatch.setattr(sfm.frames, "sample_frames", fake_sample)
+    monkeypatch.setattr(sfm, "colmap_available", lambda: True)
+
+
+def test_reconstruct_failed_no_model(monkeypatch):
+    synthetic.build(run_ffmpeg=False)
+    _fake_sampling(monkeypatch, [0, 30, 60])
+    monkeypatch.setattr(sfm, "run_colmap", lambda frames_dir, workspace: None)
+    before = _snapshot_poses()
+    result = sfm.run_reconstruct()
+    assert result["status"] == "failed_no_model"
+    assert _snapshot_poses() == before          # poses untouched
+
+
+def test_reconstruct_failed_weak(monkeypatch, tmp_path):
+    synthetic.build(run_ffmpeg=False)
+    _fake_sampling(monkeypatch, [0, 30, 60])     # 3 videos x 3 = 9 sampled
+    # model registers only one image of one video -> ~11%, < 60% and < 2 videos
+    track = synthetic.camera_track(0, 20.0, 30, 640, 360)
+    rec = dict(track[0]); rec["camera_id"] = 1; rec["name"] = "1_0.jpg"
+    cams = {1: {"width": 640, "height": 360, **synthetic.intrinsics(640, 360)}}
+    pts, cols = synthetic.generate_point_cloud()
+    model_dir = tmp_path / "weak"
+    _write_colmap_model(model_dir, [rec], cams, pts[:50], cols[:50])
+    monkeypatch.setattr(sfm, "run_colmap", lambda frames_dir, workspace: model_dir)
+    before = _snapshot_poses()
+    result = sfm.run_reconstruct()
+    assert result["status"] == "failed_weak"
+    assert _snapshot_poses() == before          # poses untouched
+
+
+def test_reconstruct_full_export_path(monkeypatch, tmp_path):
+    """The whole export path runs on injected known poses (spec M2 acceptance)."""
+    synthetic.build(run_ffmpeg=False)
+    videos = db.get_videos()
+    id_by_cam = {i: videos[i].id for i in range(len(videos))}   # cam index -> db id
+    sampled_frames = [0, 15, 30, 45, 60]
+    _fake_sampling(monkeypatch, sampled_frames)
+
+    # Build a model: videos for cam0/cam1 fully registered; cam2 registers a subset so
+    # interpolation must fill its gaps. All frame indices come from the synthetic tracks.
+    records, cameras = [], {}
+    reg_plan = {0: sampled_frames, 1: sampled_frames, 2: [0, 30, 60]}
+    for cam_idx, frame_list in reg_plan.items():
+        vid = id_by_cam[cam_idx]
+        cameras[vid] = {"width": 640, "height": 360, **synthetic.intrinsics(640, 360)}
+        track = {p["frame_idx"]: p for p in synthetic.camera_track(cam_idx, 20.0, 30, 640, 360)}
+        for f in frame_list:
+            p = track[f]
+            records.append({**p, "camera_id": vid, "name": f"{vid}_{f}.jpg"})
+    pts, cols = synthetic.generate_point_cloud()
+    model_dir = tmp_path / "good"
+    _write_colmap_model(model_dir, records, cameras, pts, cols)
+    monkeypatch.setattr(sfm, "run_colmap", lambda frames_dir, workspace: model_dir)
+
+    result = sfm.run_reconstruct()
+    assert result["status"] == "ok"
+    assert result["registered"] == 13           # 5 + 5 + 3
+    assert result["videos_in_model"] == 3
+    assert result["interpolated"] == 2           # cam2 frames 15 and 45
+
+    # cam2 poses: 5 total (3 registered + 2 interpolated), sorted, with interp flags
+    cam2_poses = db.get_poses(id_by_cam[2])
+    assert len(cam2_poses) == 5
+    assert sum(1 for p in cam2_poses if not p.registered) == 2
+    assert [p.frame_idx for p in cam2_poses] == [0, 15, 30, 45, 60]
+
+    # normalized point cloud written in the frozen PLY format; camera sphere radius ~ 10
+    assert config.POINTS_PLY.exists()
+    npoints, _ = synthetic.read_ply(config.POINTS_PLY)
+    assert len(npoints) == len(pts)
+    all_centers = []
+    for v in videos:
+        for p in db.get_poses(v.id):
+            R = quat_to_mat([p.qw, p.qx, p.qy, p.qz])
+            all_centers.append(-R.T @ np.array([p.tx, p.ty, p.tz]))
+    max_r = float(np.max(np.linalg.norm(np.array(all_centers), axis=1)))
+    assert abs(max_r - 10.0) < 0.5              # registered cams normalized to radius 10
+
+
+def test_colmap_model_parser_roundtrip(tmp_path):
+    """Writing then parsing a model recovers the poses/intrinsics/points."""
+    track = synthetic.camera_track(0, 20.0, 30, 640, 360)[:3]
+    records = [{**p, "camera_id": 1, "name": f"1_{p['frame_idx']}.jpg"} for p in track]
+    cams = {1: {"width": 640, "height": 360, **synthetic.intrinsics(640, 360)}}
+    pts, cols = synthetic.generate_point_cloud()
+    model_dir = tmp_path / "rt"
+    _write_colmap_model(model_dir, records, cams, pts[:20], cols[:20])
+
+    images = sfm.parse_images_txt(model_dir / "images.txt")
+    cameras = sfm.parse_cameras_txt(model_dir / "cameras.txt")
+    ppoints, pcolors = sfm.parse_points3d_txt(model_dir / "points3D.txt")
+    assert len(images) == 3
+    np.testing.assert_allclose(images[0]["qw"], track[0]["qw"], atol=1e-6)
+    assert cameras[1]["fx"] == pytest.approx(synthetic.intrinsics(640, 360)["fx"])
+    assert len(ppoints) == 20
+
+
+def test_parse_cameras_txt_truncated_params(tmp_path):
+    """A camera line with too few PARAMS gives a clean ValueError, not a raw IndexError."""
+    p = tmp_path / "cameras.txt"
+    p.write_text("# header\n7 PINHOLE 640 360 500.0 500.0 320.0\n")  # PINHOLE needs 4 params
+    with pytest.raises(ValueError, match="too few PARAMS"):
+        sfm.parse_cameras_txt(p)
+
+
+def test_largest_model_dir_picks_by_registered_count(tmp_path):
+    """Binary models are ranked by the images.bin header count, not file size."""
+    sparse = tmp_path / "sparse"
+    # model 0: 10 registered images but a LARGE images.bin (many keypoint observations)
+    (sparse / "0").mkdir(parents=True)
+    (sparse / "0" / "images.bin").write_bytes(struct.pack("<Q", 10) + b"\x00" * 5000)
+    # model 1: 15 registered images but a SMALL images.bin (few keypoints)
+    (sparse / "1").mkdir(parents=True)
+    (sparse / "1" / "images.bin").write_bytes(struct.pack("<Q", 15) + b"\x00" * 100)
+    assert sfm._registered_image_count(sparse / "1") == 15
+    assert sfm._largest_model_dir(sparse).name == "1"   # more images wins despite smaller file
+
+
+def test_interpolate_poses_drops_before_first_registered():
+    """Lower no-extrapolation guard: a sampled frame before the first registered one is dropped."""
+    q = [1.0, 0.0, 0.0, 0.0]
+    registered = [_pose(30, 1.0, q, [3, 0, 0]), _pose(60, 2.0, q, [6, 0, 0])]
+    sampled = [{"frame_idx": f, "t_video_s": f / 30.0} for f in (0, 30, 45, 60)]
+    out = sfm.interpolate_poses(registered, sampled)
+    frames_out = {p["frame_idx"] for p in out}
+    assert 0 not in frames_out                # tv=0 < t_first=1.0 -> dropped, no extrapolation
+    assert frames_out == {30, 45, 60}
+    by = {p["frame_idx"]: p for p in out}
+    assert by[45]["registered"] is False      # interior gap still interpolated
+
+
+def test_reconstruct_subset_leaves_others_untouched(monkeypatch, tmp_path):
+    """When only a subset of videos registers, the rest keep their existing poses (DB safety)."""
+    synthetic.build(run_ffmpeg=False)
+    videos = db.get_videos()
+    id_by_cam = {i: videos[i].id for i in range(len(videos))}
+    sampled_frames = [0, 15, 30, 45, 60]
+    _fake_sampling(monkeypatch, sampled_frames)
+
+    # Register cam0 + cam1 fully (10/15 = 67% >= 60%, 2 videos -> ok); cam2 NOT in the model.
+    records, cameras = [], {}
+    for cam_idx in (0, 1):
+        vid = id_by_cam[cam_idx]
+        cameras[vid] = {"width": 640, "height": 360, **synthetic.intrinsics(640, 360)}
+        track = {p["frame_idx"]: p for p in synthetic.camera_track(cam_idx, 20.0, 30, 640, 360)}
+        for f in sampled_frames:
+            records.append({**track[f], "camera_id": vid, "name": f"{vid}_{f}.jpg"})
+    pts, cols = synthetic.generate_point_cloud()
+    model_dir = tmp_path / "subset"
+    _write_colmap_model(model_dir, records, cameras, pts, cols)
+    monkeypatch.setattr(sfm, "run_colmap", lambda frames_dir, workspace: model_dir)
+
+    cam2_before = [(p.frame_idx, p.qw, p.tx, p.registered) for p in db.get_poses(id_by_cam[2])]
+    result = sfm.run_reconstruct()
+    assert result["status"] == "ok"
+    assert result["videos_in_model"] == 2
+    # cam2 absent from the model: its rows must be byte-identical afterward
+    cam2_after = [(p.frame_idx, p.qw, p.tx, p.registered) for p in db.get_poses(id_by_cam[2])]
+    assert cam2_after == cam2_before
+    assert len(cam2_before) == 41                       # original synthetic poses intact
+    assert len(db.get_poses(id_by_cam[0])) == 5         # cam0 replaced with sampled frames
+
+
+def test_reconstruct_no_frames_when_videos_missing(monkeypatch):
+    """All raw files absent -> videos skipped -> no_frames, poses untouched."""
+    synthetic.build(run_ffmpeg=False)
+    for f in config.RAW_DIR.glob("*.mp4"):              # clear stubs left by other tests
+        f.unlink()
+    monkeypatch.setattr(sfm, "colmap_available", lambda: True)
+    before = _snapshot_poses()
+    result = sfm.run_reconstruct()
+    assert result["status"] == "no_frames"
+    assert _snapshot_poses() == before
+
+
+def test_reconstruct_no_frames_empty_sampling(monkeypatch):
+    """Raw files exist but sampling yields nothing -> no_frames, poses untouched."""
+    synthetic.build(run_ffmpeg=False)
+    _fake_sampling(monkeypatch, [])                     # sample_frames returns []
+    before = _snapshot_poses()
+    result = sfm.run_reconstruct()
+    assert result["status"] == "no_frames"
+    assert _snapshot_poses() == before

plan/status/lane-B.md

@@ -0,0 +1,59 @@
+# Status — lane-B
+
+## Round 0 — 2026-07-16 — STATUS: ready_to_merge
+
+**Directives acknowledged:** round 0 of plan/DIRECTIVES.md (run order C→A→D→B; merge when
+acceptance passes, don't wait on siblings; maintain this status file; evidence discipline).
+
+**Environment:** COLMAP 4.1.0 (CPU-only) + ffmpeg both installed, so the pipeline was
+exercised against **real COLMAP end-to-end**, not just the synthetic export path.
+
+**Isolation note (see blocker):** worked in a dedicated git worktree
+(`/Users/m3ultra/Documents/festifun-laneB` on `lane/b-recon`) because all lane agents were
+sharing one working directory. Only the 6 lane-B-owned files were touched.
+
+**Acceptance checklist** (spec M2 + M8 geometry + lane brief) — all pass:
+- [x] `geometry.py` M8 stubs implemented: `slerp_pose` (shortest-arc, double-cover), `ray_from_pixel`
+      (COLMAP +y-down back-projection), `triangulate_rays` (closest-point + parallel guard),
+      `nearest_point_on_ray` (in-front radius cylinder). Frozen `colmap_to_threejs`/`quat_to_mat`/
+      `mat_to_quat` and the 3 frozen POSE_TEST_VECTORS untouched (`git diff 5fa7301` shows 0 removed
+      lines in test_geometry.py; the frozen funcs are unchanged).
+- [x] `frames.py`: `sharpness` (variance of Laplacian) + windowed sharpest-frame `sample_frames`.
+- [x] COLMAP TXT parsers (images/cameras/points3D) with hand-written-snippet tests + a truncated-PARAMS
+      guard (clean ValueError).
+- [x] `normalize_scene`: centroid≈0, camera-sphere radius≈10, up≈+Y — verified by `test_normalize_scene_invariants`
+      and a projection-invariance test (`test_normalize_scene_preserves_projection`).
+- [x] pose interpolation (slerp+lerp, `registered=False`, NO extrapolation past first/last registered) —
+      both guard directions tested.
+- [x] PLY export via frozen `synthetic.write_ply`; poses written via atomic `db.set_poses`.
+- [x] graceful degradation: COLMAP absent / <60% frames / <2 videos / no frames / missing files →
+      diagnostic, existing poses UNTOUCHED, DB uncorrupted. Mixed-subset case (some videos registered,
+      others kept) tested for the DB-safety invariant.
+- [x] real end-to-end: `synthetic` → `reconstruct` ran the full COLMAP CLI (feature_extractor →
+      exhaustive_matcher → mapper → image_undistorter → model_converter → my parsers), produced 6
+      components (2/15/11/16/17/37 imgs), correctly selected the largest by `images.bin` header count,
+      judged 37/120 (31%) weak, degraded gracefully, **left all 41 synthetic poses/video untouched, exit 0.**
+- [x] `pytest` green — **61 passed** (24 foundation + 37 lane-B; test_geometry.py + test_sfm.py).
+- [x] no edits outside owned files (`git diff 5fa7301 --name-only` = frames/geometry/sfm/test_geometry only;
+      untracked: test_sfm.py, this status file); no new deps.
+
+**Bugs found & fixed while validating on real COLMAP:**
+- COLMAP 4.x renamed `SiftExtraction`/`SiftMatching` → `FeatureExtraction`/`FeatureMatching`; option
+  detected from `--help` (which COLMAP prints to **stderr**) so the pipeline runs on 3.x and 4.x.
+- `_largest_model_dir` ranked binary models by file size (tracks keypoints, not image count) → now reads
+  the exact registered-image count from the `images.bin` uint64 header (validated on the real 6-component run).
+- workspace cleared before each run so a stale `database.db` can't fail reruns.
+
+**Adversarial review:** 5-lens review + verify pass (23 agents) → 18 raw findings, 5 CONFIRMED (2 code, 3
+test-coverage), all applied; 13 refuted (incl. triangulate/nearest-point "behind origin" — correct-as-written).
+
+**Blockers / questions for coordinator:**
+- **Shared working tree collision.** All four lane agents are operating in the *same* working directory
+  (`/Users/m3ultra/Documents/festifun`), not separate clones/worktrees as `plan/README.md` requires.
+  Mid-session the shared branch was switched to `lane/d-events` and the tree accumulated uncommitted edits
+  from lanes A/B/C/D at once. I isolated lane B in its own worktree so nothing is clobbered. **Recommend the
+  coordinator have each lane use `git worktree`/separate clones, and adjudicate merges centrally**, since a
+  naive `git branch -f main` from any lane could silently drop a sibling's merge. I did NOT auto-merge to
+  `main` for this reason — lane B is ready and awaiting a safe merge.
+
+**Next:** commit lane B to `lane/b-recon`; coordinator to merge to `main` (safely, given the shared-tree issue).