OpenPTV2 & Flowtracks 3D Particle Tracking Velocimetry Pipeline & AI Agent Skill

README.md

OpenPTV2 & Flowtracks 3D Particle Tracking Pipeline & AI Agent Skill

This Gist provides a complete, standalone 3D Particle Tracking Velocimetry (PTV) pipeline and AI Agent Skill using openptv2 and postptv/flowtracks.

---

📦 Contents

1. ptv_pipeline.py: A self-contained, single-file Python script to execute the complete end-to-end PTV pipeline:
   • Automated layout & calibration path verification (supports image-splitter & multi-camera setups).
   • OpenPTV2 tracking execution.
   • Flowtracks trajectory linking with exact per-trajectory velocity/acceleration differentiation.
   • Unlinked particle filtering (trajid == 0).
   • Eulerian grid binning and phase-averaged turbulence statistics (TKE, Reynolds stresses).
   • Automated value-based output validation (checks for silent-zero bugs & coordinate scaling).
2. SKILL.md: An AI Agent Skill compatible with Google Antigravity, Cursor, Claude Code, and GitHub Copilot. It equips your AI assistant with deep domain knowledge for diagnosing and orchestrating 3D PTV experiments.

---

🚀 Quick Start

1. Installation

Install the required packages (Python 3.12+):

pip install openptv2 flowtracks zarr pyyaml numpy scipy

2. Run the Pipeline

Run the pipeline on your experiment folder (containing images, calibration, and parameters_*.yaml):

python ptv_pipeline.py /path/to/experiment

Optional Arguments:

# Run specific frame range:
python ptv_pipeline.py /path/to/experiment --first 1 --last 100

# Specify frame rate (default: 1.0 fps -> velocities in mm/frame; e.g. 50.0 -> mm/s):
python ptv_pipeline.py /path/to/experiment --frame-rate 50.0

# Generate Eulerian grid with 5.0 mm voxels:
python ptv_pipeline.py /path/to/experiment --eulerian-cell 5.0

# Target a specific parameters file:
python ptv_pipeline.py /path/to/experiment --params parameters_Run1.yaml

---

🤖 Using the AI Agent Skill (SKILL.md)

To equip your AI assistant (e.g. Antigravity, Cursor, Claude Code) with 3D PTV domain knowledge:

• Google Antigravity / Claude Code / Agent Workspaces: Save SKILL.md into your workspace under:

  .agents/skills/openptv/SKILL.md
  

• Your assistant will automatically use the skill whenever you ask to run PTV processing, diagnose tracking failures, verify calibration, or post-process velocity fields.

ptv_pipeline.py

#!/usr/bin/env python3
"""Standalone 3D Particle Tracking Velocimetry (PTV) Pipeline.

Uses openptv2 and postptv (flowtracks) to process raw PTV experiment folders
into 3D trajectories, velocity/acceleration fields, and Eulerian statistics.

Usage:
    python ptv_pipeline.py <experiment_dir> [options]

Options:
    --first <int>           First frame to process
    --last <int>            Last frame to process
    --params <filename>     Target parameters YAML file
    --frame-rate <float>    Camera acquisition frame rate in Hz (default: 1.0)
    --eulerian-cell <float> Eulerian voxel cell size in mm (optional)
    --min-length <int>      Minimum trajectory frame length (default: 3)
"""

from __future__ import annotations

import argparse
import sys
from pathlib import Path
from typing import Any

import numpy as np
import yaml

INPUT_PATH_KEYS = [
    ("cal_ori", "img_ori", True),
    ("cal_ori", "img_cal_name", True),
    ("cal_ori", "fixp_name", False),
]


def is_placeholder(val: str) -> bool:
    """Check if a calibration string is an image-splitter placeholder like '---'."""
    s = val.strip()
    return not s or set(s) <= {"-"}


def find_params_file(data_dir: Path, requested_name: str | None = None) -> Path:
    """Locate the parameters_*.yaml file in the experiment directory."""
    if requested_name:
        p = data_dir / requested_name
        if p.is_file():
            return p
        matches = list(data_dir.glob(requested_name))
        if matches:
            return matches[0]
        raise FileNotFoundError(
            f"Parameters file '{requested_name}' not found in {data_dir}"
        )

    candidates = sorted(data_dir.glob("parameters_*.yaml"))
    if not candidates:
        candidates = sorted(data_dir.glob("parameters*.yaml"))
    if not candidates:
        raise FileNotFoundError(f"No parameters_*.yaml found in {data_dir}")

    # Prefer batch over sample if multiple exist
    batches = [c for c in candidates if "_batch" in c.stem]
    return batches[0] if batches else candidates[0]


def check_inputs(
    data_dir: Path, params_file: Path, first: int, last: int
) -> dict[str, Any]:
    """Verify that all calibration files and frame images referenced in YAML exist."""
    params = yaml.safe_load(params_file.read_text(encoding="utf-8")) or {}
    print(f"[pipeline] verifying inputs from {params_file.name}...")

    missing = []
    for section, key, is_list in INPUT_PATH_KEYS:
        val = (params.get(section) or {}).get(key)
        if val is None:
            continue
        items = val if is_list else [val]
        for item in items:
            if isinstance(item, str) and item.strip() and not is_placeholder(item):
                resolved = (data_dir / item).resolve()
                if not resolved.exists():
                    missing.append(f"  MISSING: {item} -> {resolved}")
                else:
                    print(f"  [ok] {item}")

    # Verify endpoint images
    base_names = ((params.get("sequence") or {}).get("base_name")) or []
    for pattern in base_names:
        if (
            not isinstance(pattern, str)
            or "%" not in pattern
            or is_placeholder(pattern)
        ):
            continue
        for frame in {first, last}:
            rel = pattern % frame
            resolved = (data_dir / rel).resolve()
            if not resolved.exists():
                missing.append(f"  MISSING image: {rel} -> {resolved}")
            else:
                print(f"  [ok] image frame {frame}: {rel}")

    if missing:
        raise FileNotFoundError("Missing required input files:\n" + "\n".join(missing))
    return params


def compute_derivatives(
    pos: np.ndarray, time: np.ndarray, dt: float
) -> tuple[np.ndarray, np.ndarray]:
    """Compute central difference velocity and acceleration for a single trajectory."""
    n = len(pos)
    vel = np.zeros_like(pos, dtype=np.float64)
    accel = np.zeros_like(pos, dtype=np.float64)

    if n < 2:
        return vel, accel

    if n == 2:
        d_time = (time[1] - time[0]) * dt
        if d_time > 0:
            vel[:] = (pos[1] - pos[0]) / d_time
        return vel, accel

    # First point (forward difference)
    dt0 = (time[1] - time[0]) * dt
    if dt0 > 0:
        vel[0] = (pos[1] - pos[0]) / dt0

    # Interior points (central difference)
    dt_mid = (time[2:] - time[:-2]) * dt
    valid_mid = dt_mid > 0
    if np.any(valid_mid):
        vel[1:-1][valid_mid] = (
            pos[2:][valid_mid] - pos[:-2][valid_mid]
        ) / dt_mid[valid_mid, None]

    # Last point (backward difference)
    dtn = (time[-1] - time[-2]) * dt
    if dtn > 0:
        vel[-1] = (pos[-1] - pos[-2]) / dtn

    # Acceleration from velocities (central difference)
    if n >= 3:
        if dt0 > 0:
            accel[0] = (vel[1] - vel[0]) / dt0
        accel[1:-1][valid_mid] = (
            vel[2:][valid_mid] - vel[:-2][valid_mid]
        ) / dt_mid[valid_mid, None]
        if dtn > 0:
            accel[-1] = (vel[-1] - vel[-2]) / dtn

    return vel, accel


def run_openptv_tracking(
    data_dir: Path, params_file: Path, first: int, last: int
) -> None:
    """Execute OpenPTV2 tracking on the experiment folder."""
    try:
        from openptv2.main import main as optv_main

        print(
            f"[pipeline] running openptv2 tracking on {data_dir.name} "
            f"(frames {first}..{last})..."
        )
        optv_main([str(data_dir), str(params_file.name), str(first), str(last)])
    except ImportError:
        import subprocess

        print("[pipeline] calling openptv2 via subprocess...")
        res = subprocess.run(
            [
                sys.executable,
                "-m",
                "openptv2",
                str(data_dir),
                str(params_file.name),
                str(first),
                str(last),
            ],
            cwd=str(data_dir),
            capture_output=True,
            text=True,
        )
        if res.returncode != 0:
            raise RuntimeError(f"openptv2 failed:\n{res.stderr}\n{res.stdout}")


def build_trajectories_table(
    res_dir: Path, frame_rate: float = 1.0
) -> dict[str, np.ndarray]:
    """Read ptv_is linkage from res/ and build structured trajectory arrays."""
    from flowtracks.io import trajectories_ptvis

    ptv_files = sorted(res_dir.glob("ptv_is.*"))
    if not ptv_files:
        raise FileNotFoundError(f"No ptv_is.* output files found in {res_dir}")

    first_frame = int(ptv_files[0].suffix.lstrip("."))
    last_frame = int(ptv_files[-1].suffix.lstrip("."))
    print(
        f"[pipeline] reading trajectories from ptv_is files "
        f"({first_frame}..{last_frame})..."
    )

    # Read trajectories using flowtracks
    trajects = trajectories_ptvis(str(res_dir), first=first_frame, last=last_frame)
    dt = 1.0 / frame_rate if frame_rate > 0 else 1.0

    all_pos, all_vel, all_accel, all_time, all_trajid = [], [], [], [], []

    for tr in trajects:
        tr_id = int(tr.trajid())
        pos = np.asarray(tr.pos(), dtype=np.float64)
        time = np.asarray(tr.time(), dtype=np.int64)

        if tr_id == 0 or len(pos) == 0:
            # Unlinked particles: no valid derivatives
            vel = np.zeros_like(pos)
            accel = np.zeros_like(pos)
        else:
            sort_idx = np.argsort(time)
            pos = pos[sort_idx]
            time = time[sort_idx]
            vel, accel = compute_derivatives(pos, time, dt)

        all_pos.append(pos)
        all_vel.append(vel)
        all_accel.append(accel)
        all_time.append(time)
        all_trajid.append(np.full(len(pos), tr_id, dtype=np.int64))

    if not all_pos:
        raise ValueError("No particle data was recovered from tracking output.")

    pos_cat = np.concatenate(all_pos, axis=0)
    vel_cat = np.concatenate(all_vel, axis=0)
    accel_cat = np.concatenate(all_accel, axis=0)
    time_cat = np.concatenate(all_time, axis=0)
    trajid_cat = np.concatenate(all_trajid, axis=0)

    return {
        "pos": pos_cat,
        "vel": vel_cat,
        "accel": accel_cat,
        "time": time_cat,
        "trajid": trajid_cat,
    }


def save_zarr(res_dir: Path, data: dict[str, np.ndarray]) -> Path:
    """Save structured trajectory data into a Zarr store."""
    import zarr

    zarr_path = res_dir / "run.zarr"
    store = zarr.open_group(str(zarr_path), mode="a")
    traj_grp = store.require_group("trajectories")

    for key, arr in data.items():
        if key in traj_grp:
            del traj_grp[key]
        traj_grp.create_dataset(key, data=arr, overwrite=True)

    print(f"[pipeline] saved structured trajectories to {zarr_path}/trajectories")
    return zarr_path


def compute_eulerian_grid(
    pos: np.ndarray,
    vel: np.ndarray,
    cell_size: float,
    min_samples: int = 3,
) -> dict[str, Any]:
    """Bin 3D velocities onto a regular Eulerian Cartesian grid."""
    # Use percentiles to avoid outlier explosion
    p_low = np.percentile(pos, 0.5, axis=0)
    p_high = np.percentile(pos, 99.5, axis=0)

    x_edges = np.arange(p_low[0], p_high[0] + cell_size, cell_size)
    y_edges = np.arange(p_low[1], p_high[1] + cell_size, cell_size)
    z_edges = np.arange(p_low[2], p_high[2] + cell_size, cell_size)

    nx, ny, nz = len(x_edges) - 1, len(y_edges) - 1, len(z_edges) - 1
    counts = np.zeros((nx, ny, nz), dtype=np.int32)
    u_sum = np.zeros((nx, ny, nz, 3), dtype=np.float64)
    u_sq_sum = np.zeros((nx, ny, nz, 3), dtype=np.float64)

    ix = np.clip(np.digitize(pos[:, 0], x_edges) - 1, 0, nx - 1)
    iy = np.clip(np.digitize(pos[:, 1], y_edges) - 1, 0, ny - 1)
    iz = np.clip(np.digitize(pos[:, 2], z_edges) - 1, 0, nz - 1)

    for i in range(len(pos)):
        cx, cy, cz = ix[i], iy[i], iz[i]
        counts[cx, cy, cz] += 1
        v = vel[i]
        u_sum[cx, cy, cz] += v
        u_sq_sum[cx, cy, cz] += v * v

    valid = counts >= min_samples
    u_mean = np.full((nx, ny, nz, 3), np.nan)
    u_mean[valid] = u_sum[valid] / counts[valid, None]

    # Reynolds stresses: <u'u'> = <u^2> - <u>^2
    reynolds_diag = np.full((nx, ny, nz, 3), np.nan)
    reynolds_diag[valid] = (u_sq_sum[valid] / counts[valid, None]) - (
        u_mean[valid] ** 2
    )
    reynolds_diag[reynolds_diag < 0] = 0.0

    # TKE = 0.5 * (u'^2 + v'^2 + w'^2)
    tke = np.full((nx, ny, nz), np.nan)
    tke[valid] = 0.5 * np.sum(reynolds_diag[valid], axis=1)

    occupancy = float(np.sum(valid) / counts.size * 100.0)
    print(
        f"[eulerian] grid ({nx}, {ny}, {nz}) -> "
        f"{occupancy:.1f}% occupied voxels (>= {min_samples} samples)"
    )

    return {
        "x": 0.5 * (x_edges[:-1] + x_edges[1:]),
        "y": 0.5 * (y_edges[:-1] + y_edges[1:]),
        "z": 0.5 * (z_edges[:-1] + z_edges[1:]),
        "counts": counts,
        "mean_velocity": u_mean,
        "reynolds_diag": reynolds_diag,
        "tke": tke,
    }


def verify_results(data: dict[str, np.ndarray]) -> None:
    """THE FIRST RULE: Assert on actual numerical values produced by the pipeline."""
    pos = data["pos"]
    vel = data["vel"]
    trajid = data["trajid"]

    if len(pos) == 0:
        raise ValueError("Pipeline produced 0 total particle points!")

    n_trajs = len(np.unique(trajid[trajid > 0]))
    n_unlinked = int((trajid == 0).sum())
    speeds = np.linalg.norm(vel[trajid > 0], axis=1) if n_trajs > 0 else np.array([])
    valid_speeds = speeds[speeds > 0]

    print("\n" + "=" * 50)
    print("           PIPELINE QUALITY AUDIT           ")
    print("=" * 50)
    print(f"Total points       : {len(pos):,}")
    print(f"Linked trajectories: {n_trajs:,}")
    print(
        f"Unlinked fraction  : {n_unlinked / len(pos) * 100:.1f}% "
        f"({n_unlinked:,} rows)"
    )

    if len(valid_speeds) == 0:
        raise ValueError("CRITICAL BUG: All trajectory velocities are exactly 0.0!")

    print(f"Mean speed (|v|)   : {valid_speeds.mean():.3f}")
    print(f"Max speed (|v|)    : {valid_speeds.max():.3f}")
    print(f"Coordinate extent  : X=[{pos[:, 0].min():.1f}, {pos[:, 0].max():.1f}]")
    print(f"                     Y=[{pos[:, 1].min():.1f}, {pos[:, 1].max():.1f}]")
    print(f"                     Z=[{pos[:, 2].min():.1f}, {pos[:, 2].max():.1f}]")
    print(f"Finite values check: {np.isfinite(vel).all()}")
    print("=" * 50 + "\n")


def main() -> None:
    parser = argparse.ArgumentParser(description="Standalone 3D OpenPTV Pipeline")
    parser.add_argument("path", type=Path, help="Experiment directory path")
    parser.add_argument("--first", type=int, help="First frame number")
    parser.add_argument("--last", type=int, help="Last frame number")
    parser.add_argument("--params", type=str, help="Parameters YAML filename")
    parser.add_argument(
        "--frame-rate",
        type=float,
        default=1.0,
        help="Camera acquisition frame rate in Hz",
    )
    parser.add_argument(
        "--eulerian-cell", type=float, help="Eulerian cell voxel size in mm"
    )
    args = parser.parse_args()

    data_dir = args.path.resolve()
    if not data_dir.is_dir():
        raise SystemExit(f"Directory not found: {data_dir}")

    params_file = find_params_file(data_dir, args.params)
    params = yaml.safe_load(params_file.read_text(encoding="utf-8")) or {}

    first = (
        args.first
        if args.first is not None
        else int((params.get("sequence") or {}).get("first", 1))
    )
    last = (
        args.last
        if args.last is not None
        else int((params.get("sequence") or {}).get("last", 10))
    )

    # 1. Preflight check
    check_inputs(data_dir, params_file, first, last)

    # 2. Tracking execution
    run_openptv_tracking(data_dir, params_file, first, last)

    # 3. Post-process trajectories and compute derivatives
    res_dir = data_dir / "res"
    traj_data = build_trajectories_table(res_dir, frame_rate=args.frame_rate)

    # 4. Save to Zarr
    save_zarr(res_dir, traj_data)

    # 5. Eulerian Grid (optional)
    if args.eulerian_cell:
        linked_mask = traj_data["trajid"] > 0
        grid_data = compute_eulerian_grid(
            traj_data["pos"][linked_mask],
            traj_data["vel"][linked_mask],
            cell_size=args.eulerian_cell,
        )
        import zarr

        store = zarr.open_group(str(res_dir / "run.zarr"), mode="a")
        eul_grp = store.require_group("eulerian")
        for k, v in grid_data.items():
            if k in eul_grp:
                del eul_grp[k]
            eul_grp.create_dataset(k, data=v, overwrite=True)
        print(f"[pipeline] saved Eulerian grid to {res_dir / 'run.zarr'}/eulerian")

    # 6. Quality self-check
    verify_results(traj_data)
    print("[pipeline] SUCCESS: 3D PTV pipeline completed cleanly.")


if __name__ == "__main__":
    main()

SKILL.md

name	openptv
description	Run, process, and troubleshoot 3D Particle Tracking Velocimetry (PTV) experiments with openptv2 and postptv (flowtracks). Use when processing calibration, running particle tracking, building 3D trajectories, computing velocity/acceleration derivatives, generating Eulerian grids, phase-averaging, or diagnosing empty/silent errors.

openptv: 3D Particle Tracking Velocimetry (PTV)

Turns an experiment directory (images + calibration + an openptv2 parameters_*.yaml) into 3D particle trajectories, velocity/acceleration fields, and Eulerian turbulence statistics using openptv2 and postptv (flowtracks).

---

1. Prerequisites

pip install openptv2 flowtracks zarr pyyaml numpy scipy

---

2. Core Domain Knowledge

A. openptv2 Configuration (parameters_*.yaml)

• Relative Calibration Paths: parameters_*.yaml paths resolve relative to the folder being processed. In multi-run setups, subfolder YAMLs reach shared calibration as ../calibration/.... Never flatten this folder structure.
• Image-Splitter Mode (splitter: true or cal_splitter: true):
  • 4 "cameras" are quadrants of one image file.
  • sequence.base_name repeats the same file path four times.
  • Unused cal_ori.img_cal_name slots hold '---' placeholders. These dashes are not filenames — do not treat them as missing files.
• Coordinates & Units: openptv2 writes in millimetres (mm) by default. Always verify coordinate extents before sizing Eulerian grids.
• Output Artifacts: openptv2 generates *_targets and ptv_is.* files during tracking.

B. postptv & flowtracks (Trajectories & Derivatives)

• Catch-all Bucket (trajid == 0): trajid == 0 is the collection of unlinked particles, NOT a trajectory. It contains multiple particles sharing the same frames. Differentiating it causes $\Delta t = 0$ division errors ($NaN$ / $\infty$). Always filter trajid > 0 for velocity/acceleration statistics.
• Numerical Derivatives: Trajectories may contain frame gaps. Derivatives must be calculated per trajectory along its sorted frame sequence using central difference: $$v_i = \frac{x_{i+1} - x_{i-1}}{t_{i+1} - t_{i-1}}$$
• Storage: Structured trajectories are stored in Zarr format (res/run.zarr/trajectories) with flat 1D/2D arrays: pos $(N, 3)$, vel $(N, 3)$, accel $(N, 3)$, time $(N,)$, and trajid $(N,)$.

---

3. THE FIRST RULE: A "Successful" Run Can Still Be Empty

OpenPTV pipelines can execute without raising exceptions while producing zero or invalid data. Always assert on the values:

import numpy as np, zarr

g = zarr.open_group("res/run.zarr", mode="r")["trajectories"]
pos = np.asarray(g["pos"][:])
vel = np.asarray(g["vel"][:])
trajid = np.asarray(g["trajid"][:])
speed = np.linalg.norm(vel, axis=1)

print(f"Trajectories: {len(np.unique(trajid[trajid > 0]))}, Total points: {len(pos)}")
print(f"Valid velocity rows: {int((speed > 0).sum())}, Unlinked: {int((trajid == 0).sum())}")
print(f"Mean speed: {speed[speed > 0].mean():.3f}, Finite check: {np.isfinite(vel).all()}")
print(f"Position bounds: {pos.min(0).round(2)} to {pos.max(0).round(2)}")

Diagnostic Table

Symptom	Probable Cause	Corrective Action
Every velocity is exactly 0	Derivatives were not computed or unlinked rows included	Calculate central differences per trajectory (trajid > 0).
inf / NaN velocities	trajid == 0 was differentiated across duplicate frame numbers	Filter out trajid == 0 before differentiation.
Eulerian field all zeros	Grid bounding box does not overlap particles	Check coordinate units (mm vs m); compute bounds from particle percentiles.
Position extent ~0.05	Coordinates in meters, but grid assumes millimeters	Inspect calibration units; convert units consistently.
High unlinked fraction (>50%)	Search radius / intensity threshold misconfigured	Adjust track.yaml (search radius, acceleration limit).

---

4. Eulerian Grid Sizing & Post-Processing

Never guess the Eulerian grid extent. Estimate it from particle positions:

1. Compute percentiles: Use p0.5 to p99.5 rather than min/max to exclude stray outliers.
2. Select cell size:
   • $5\text{ mm}$ cells $\rightarrow$ higher sample count per voxel ($\sim 30\text{--}60%$ occupancy).
   • $2\text{ mm}$ cells $\rightarrow$ finer spatial resolution but sparser occupancy.
3. Compute Turbulence Statistics:
   • Mean Velocity: $\bar{\mathbf{u}} = \frac{1}{N} \sum \mathbf{u}_i$
   • Reynolds Stresses: $\overline{u'_i u'_j} = \frac{1}{N} \sum (u_i - \bar{u}_i)(u_j - \bar{u}_j)$
   • Turbulent Kinetic Energy (TKE): $k = \frac{1}{2}(\overline{u'^2} + \overline{v'^2} + \overline{w'^2})$

---

5. End-to-End Execution Script

Use ptv_pipeline.py to run the full automated workflow on any experiment folder:

# Standard run
python ptv_pipeline.py /path/to/experiment

# With 50 FPS and 5mm Eulerian grid
python ptv_pipeline.py /path/to/experiment --frame-rate 50.0 --eulerian-cell 5.0