From 6487919bd1518fe05a96959235e286e0eaff5028 Mon Sep 17 00:00:00 2001 From: Alexander Refsum Jensenius Date: Thu, 16 Jul 2026 09:43:39 +0200 Subject: [PATCH 1/3] Add _posetools: landmark-trajectory pose extraction and derived motion signals New module implementing the array-level pose workflow used across the author's sound-motion studies: video -> tidy per-landmark trajectory arrays/CSV -> derived motion signals. Complements the rendering-oriented MgVideo.pose() pipeline and the per-frame PoseEstimator interface. - extract_pose_landmarks: consolidated MediaPipe extractor (FFmpeg-piped frame reader, fps resampling, resize, model complexity, optional 3D world landmarks, confidence thresholds, NaN on detection failure, detection-rate summary, optional tidy CSV). MediaPipe is imported lazily (optional [pose] extra, already declared in pyproject); both the legacy Solutions API (0.10.14-era wheels) and the Tasks API PoseLandmarker in VIDEO mode (newer 0.10.x wheels) are supported, sharing the model cache with MgVideo.pose(). - midpoint: NaN-propagating trajectory midpoint (e.g. shoulder midpoint as upper-torso proxy). - limb_speed_from_landmarks: confidence-gated central-difference pixel speed with bilateral (element-wise max) limb merge and NaN-aware smoothing; documents the speed-peak-precedes-contact timing bias. - impact_events: candidate impact detection from point trajectories via double-central-difference acceleration magnitude, bilateral max, and an inline relative-threshold peak picker (provisional defaults from the cymbal study; the general picker lives in the sibling _peaks PR). Consolidated/ported from the stillstanding (mp_extract_westney.py, pose_motion.py), Westney-comparisons (concert_mediapipe.py, reh_pose.py, a1_labstage.py) and cymbal-comparison study code/papers (all author-owned analysis code). Multi-view Umeyama pose fusion is deliberately deferred to a phase-2 follow-up. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_0147io2kMk8M6jcNFNt1G96m --- musicalgestures/__init__.py | 9 + musicalgestures/_posetools.py | 563 ++++++++++++++++++++++++++++++++++ 2 files changed, 572 insertions(+) create mode 100644 musicalgestures/_posetools.py diff --git a/musicalgestures/__init__.py b/musicalgestures/__init__.py index 03e4fe8..ac2ab64 100644 --- a/musicalgestures/__init__.py +++ b/musicalgestures/__init__.py @@ -69,3 +69,12 @@ def __init__(self): rayleigh_test, synchrony, ) +# Landmark-trajectory pose tools. extract_pose_landmarks needs the optional +# mediapipe package ([pose] extra) but imports it lazily, so this is safe +# without it; the other helpers are numpy-only. +from musicalgestures._posetools import ( + extract_pose_landmarks, + midpoint, + limb_speed_from_landmarks, + impact_events, +) diff --git a/musicalgestures/_posetools.py b/musicalgestures/_posetools.py new file mode 100644 index 0000000..108d711 --- /dev/null +++ b/musicalgestures/_posetools.py @@ -0,0 +1,563 @@ +""" +Landmark-trajectory pose tools. + +This module implements the *array-level* pose workflow used in several of the +fourMs sound--motion studies: video file -> tidy per-landmark trajectory +arrays (and optionally CSV) -> derived motion signals (limb speed, impact +events). + +It complements — and does not replace — the rendering-oriented +``MgVideo.pose()`` pipeline in :mod:`musicalgestures._pose` (overlaid skeleton +video, average-pose image, trajectory image, keypoint CSV) and the per-frame +:class:`musicalgestures._pose_estimator.PoseEstimator` interface. Use this +module when you want plain numpy trajectories for downstream signal analysis +(quantity of motion, cross-modal alignment, event detection) rather than +rendered output. + +Only :func:`extract_pose_landmarks` needs MediaPipe (an optional dependency, +imported lazily). The derived-signal helpers (:func:`midpoint`, +:func:`limb_speed_from_landmarks`, :func:`impact_events`) are numpy-only and +also work on landmark/point trajectories from any other source (OpenPose, +YOLO-pose, motion capture). + +Landmark indices follow the 33-landmark MediaPipe Pose (BlazePose GHUM) +topology used by the mediapipe 0.10.x wheels (e.g. 0 = nose, 11/12 = +left/right shoulder, 13/14 = elbows, 15/16 = wrists); see +:data:`musicalgestures._pose_estimator.MEDIAPIPE_LANDMARK_NAMES` for the full +index -> name mapping. +""" + +from __future__ import annotations + +import os +import subprocess +import warnings + +import numpy as np + +from musicalgestures._utils import get_widthheight, get_fps + +# Landmark names are defined next to the per-frame estimator so both pose +# workflows agree on the index -> name mapping. _pose_estimator imports +# mediapipe lazily, so this import is safe without mediapipe installed. +from musicalgestures._pose_estimator import MEDIAPIPE_LANDMARK_NAMES + + +def extract_pose_landmarks( + filename: str, + fps: float | None = None, + width: int | None = None, + model_complexity: int = 1, + world_landmarks: bool = False, + min_detection_confidence: float = 0.5, + min_tracking_confidence: float = 0.5, + max_frames: int | None = None, + target_name: str | None = None, + quiet: bool = True, + verbose: bool = True) -> dict: + """ + Run MediaPipe Pose over a whole video and return tidy per-landmark trajectories. + + The video is decoded through an FFmpeg raw-video pipe (optionally resampled + to a lower frame rate and resized), each frame is passed to MediaPipe Pose, + and the 33 landmarks are collected into plain numpy arrays: pixel + coordinates in the analysis frame plus the per-landmark ``visibility`` + score, with all-NaN rows on frames where no pose was detected, and a + detection-rate summary. This is the consolidated version of several + near-identical study extractors; downstream method choices (filtering, QoM, + alignment) are deliberately *not* baked in here. + + MediaPipe is an optional dependency (``pip install musicalgestures[pose]``) + and is imported lazily, so importing this module works without it. Both + mediapipe API families are supported: the legacy Solutions API + (``mp.solutions.pose.Pose``, wheels up to ~0.10.14) and the Tasks API + (``PoseLandmarker`` in VIDEO running mode, newer 0.10.x wheels where the + Solutions API was removed). With the Tasks API the model file is + auto-downloaded and cached in ``musicalgestures/models/`` on first use + (shared with ``MgVideo.pose()``). + + Args: + filename (str): Path to the input video file. + fps (float, optional): Analysis frame rate. Frames are resampled to this + rate by FFmpeg before pose estimation (e.g. 12.5 to halve a 25 fps + video). Defaults to None (native frame rate). + width (int, optional): Resize the analysis frames to this width in + pixels, keeping the aspect ratio. Smaller frames are much faster and + are usually sufficient for trajectory-level analysis (the studies + used 256-640 px). Defaults to None (native resolution). + model_complexity (int, optional): MediaPipe model variant: 0 (lite), + 1 (full) or 2 (heavy). Defaults to 1. + world_landmarks (bool, optional): Whether to also collect MediaPipe's 3D + world landmarks (metres, hip-centred). Defaults to False. + min_detection_confidence (float, optional): MediaPipe person-detection + confidence threshold (also used as the presence threshold with the + Tasks API). Defaults to 0.5. + min_tracking_confidence (float, optional): MediaPipe landmark-tracking + confidence threshold. Defaults to 0.5. + max_frames (int, optional): Stop after this many analysed frames (handy + for quick tests). Defaults to None (whole video). + target_name (str, optional): If given, also write the trajectories to + this path as a tidy CSV with columns ``time`` and, per landmark + name, ``_x``, ``_y``, ``_v`` (and ``_wx``, + ``_wy``, ``_wz`` when ``world_landmarks=True``). + Defaults to None (no file written). + quiet (bool, optional): Suppress MediaPipe's native C++/GL console logs + during inference. Defaults to True. + verbose (bool, optional): Print a one-line detection-rate summary per + video. Defaults to True. + + Returns: + dict: A dictionary with keys: + + - ``time`` (np.ndarray, shape (F,)): Frame timestamps in seconds. + - ``landmarks`` (np.ndarray, shape (F, 33, 3)): Per frame and + landmark ``(x_px, y_px, visibility)`` in analysis-frame pixels; + all-NaN rows where no pose was detected. + - ``world`` (np.ndarray, shape (F, 33, 3) or None): 3D world + landmarks ``(x, y, z)`` in metres (hip-centred) when + ``world_landmarks=True``, else None. + - ``detected`` (np.ndarray of bool, shape (F,)): Per-frame + detection flags. + - ``detection_rate`` (float): Fraction of frames with a detected + pose. + - ``fps`` (float): Analysis frame rate of the returned arrays. + - ``width``, ``height`` (int): Analysis frame size in pixels. + - ``names`` (list of str): The 33 landmark names (row order of the + landmark axis). + + Source: + Consolidated from the author's study extractors: stillstanding + (mp_extract_westney.py, pose_motion.py) and Westney-comparisons + (concert_mediapipe.py, reh_pose.py, a1_labstage.py) (Jensenius). + """ + try: + import mediapipe as mp + except ImportError as exc: + raise ImportError( + "MediaPipe is required for extract_pose_landmarks() but is not installed. " + "Install the optional pose dependencies with: pip install musicalgestures[pose]" + ) from exc + + # fd-level stderr silencer shared with MgVideo.pose() (MediaPipe's C++/GL + # logs bypass Python logging). + from musicalgestures._pose import _suppress_native_stderr + + n_landmarks = len(MEDIAPIPE_LANDMARK_NAMES) + + # --- probe the source and build the FFmpeg decode pipe ------------------- + w0, h0 = get_widthheight(filename) + native_fps = get_fps(filename) + sample_fps = float(fps) if fps else float(native_fps) + if width: + w = int(width) + h = int(round(h0 * w / w0)) + else: + w, h = int(w0), int(h0) + + vf = [] + if fps: + vf.append(f"fps={sample_fps}") + if width: + vf.append(f"scale={w}:{h}") + cmd = ["ffmpeg", "-v", "error", "-i", filename] + if vf: + cmd += ["-vf", ",".join(vf)] + cmd += ["-pix_fmt", "rgb24", "-f", "rawvideo", "-"] + + # --- backend: legacy Solutions API or Tasks API --------------------------- + # The study scripts were written against the legacy Solutions API + # (mediapipe 0.10.14); newer 0.10.x wheels (e.g. 0.10.35, as used in the + # cymbal study) removed it in favour of the Tasks API. Support both. + use_solutions = hasattr(mp, "solutions") + + def _read_result(lms, wlms): + """Convert one frame's MediaPipe landmark lists to (33, 3) arrays.""" + a = np.array([[l.x * w, l.y * h, l.visibility] for l in lms], dtype=np.float64) + if world_landmarks: + b = (np.array([[l.x, l.y, l.z] for l in wlms], dtype=np.float64) + if wlms is not None else np.full((n_landmarks, 3), np.nan)) + else: + b = None + return a, b + + times, lm2d, lm3d, detected = [], [], [], [] + frame_bytes = w * h * 3 + stopped_early = False + proc = subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE) + try: + with _suppress_native_stderr(quiet): + if use_solutions: + pose_ctx = mp.solutions.pose.Pose( + static_image_mode=False, + model_complexity=model_complexity, + min_detection_confidence=min_detection_confidence, + min_tracking_confidence=min_tracking_confidence) + close_backend = pose_ctx.close + + def _process(frame_rgb, t): + res = pose_ctx.process(frame_rgb) + if res.pose_landmarks: + return _read_result( + res.pose_landmarks.landmark, + res.pose_world_landmarks.landmark + if getattr(res, "pose_world_landmarks", None) else None) + return None, None + else: + # Tasks API: reuse the model download/cache logic of the + # per-frame estimator so both pose workflows share one model + # file in musicalgestures/models/. + from musicalgestures._pose_estimator import MediaPipePoseEstimator + model_path = MediaPipePoseEstimator( + model_complexity=model_complexity)._get_model_path() + BaseOptions = mp.tasks.BaseOptions + options = mp.tasks.vision.PoseLandmarkerOptions( + base_options=BaseOptions(model_asset_path=str(model_path)), + running_mode=mp.tasks.vision.RunningMode.VIDEO, + min_pose_detection_confidence=min_detection_confidence, + min_pose_presence_confidence=min_detection_confidence, + min_tracking_confidence=min_tracking_confidence) + landmarker = mp.tasks.vision.PoseLandmarker.create_from_options(options) + close_backend = landmarker.close + + def _process(frame_rgb, t): + img = mp.Image(image_format=mp.ImageFormat.SRGB, data=frame_rgb) + res = landmarker.detect_for_video(img, int(round(t * 1000))) + if res.pose_landmarks: + return _read_result( + res.pose_landmarks[0], + res.pose_world_landmarks[0] + if res.pose_world_landmarks else None) + return None, None + + try: + fi = 0 + while True: + if max_frames is not None and fi >= max_frames: + stopped_early = True + proc.terminate() + break + buf = proc.stdout.read(frame_bytes) + if buf is None or len(buf) < frame_bytes: + break + frame = np.frombuffer(buf, np.uint8).reshape(h, w, 3) + t = fi / sample_fps + a, b = _process(frame, t) + if a is None: + a = np.full((n_landmarks, 3), np.nan) + b = np.full((n_landmarks, 3), np.nan) if world_landmarks else None + detected.append(False) + else: + detected.append(True) + times.append(t) + lm2d.append(a) + if world_landmarks: + lm3d.append(b) + fi += 1 + finally: + close_backend() + finally: + proc.stdout.close() + err = proc.stderr.read().decode(errors="replace").strip() + proc.stderr.close() + proc.wait() + # Stopping at max_frames breaks FFmpeg's output pipe on purpose, so + # suppress its resulting broken-pipe complaints. + if err and not stopped_early: + print(f"FFmpeg warnings while decoding {filename}:\n{err}") + + n_frames = len(times) + result = { + "time": np.asarray(times, dtype=np.float64), + "landmarks": (np.asarray(lm2d, dtype=np.float64) + if n_frames else np.empty((0, n_landmarks, 3))), + "world": ((np.asarray(lm3d, dtype=np.float64) + if n_frames else np.empty((0, n_landmarks, 3))) + if world_landmarks else None), + "detected": np.asarray(detected, dtype=bool), + "detection_rate": float(np.mean(detected)) if n_frames else 0.0, + "fps": sample_fps, + "width": w, + "height": h, + "names": list(MEDIAPIPE_LANDMARK_NAMES), + } + + if verbose: + print(f"{os.path.basename(filename)}: {n_frames} frames at " + f"{sample_fps:g} fps ({w}x{h}), pose detected in " + f"{100.0 * result['detection_rate']:.0f}% of frames.") + + if target_name is not None: + _write_landmarks_csv(target_name, result) + + return result + + +def _write_landmarks_csv(path: str, result: dict) -> None: + """Write an extract_pose_landmarks() result dict as a tidy CSV file.""" + names = result["names"] + cols = [result["time"][:, None]] + header = ["time"] + lm = result["landmarks"] + for j, name in enumerate(names): + cols.append(lm[:, j, :]) + header += [f"{name}_x", f"{name}_y", f"{name}_v"] + if result["world"] is not None: + wl = result["world"] + for j, name in enumerate(names): + cols.append(wl[:, j, :]) + header += [f"{name}_wx", f"{name}_wy", f"{name}_wz"] + data = np.hstack(cols) if len(result["time"]) else np.empty((0, len(header))) + np.savetxt(path, data, delimiter=",", header=",".join(header), comments="") + + +def midpoint(a: np.ndarray, b: np.ndarray) -> np.ndarray: + """ + Element-wise midpoint of two landmark trajectories. + + Typical use is the shoulder midpoint (MediaPipe landmarks 11 and 12) as an + upper-torso proxy point, e.g. ``midpoint(lm[:, 11, :2], lm[:, 12, :2])``. + NaNs (detection dropouts) propagate: the midpoint is NaN wherever either + input is NaN. + + Args: + a (np.ndarray): First trajectory, any shape (e.g. (F, 2)). + b (np.ndarray): Second trajectory, broadcast-compatible with ``a``. + + Returns: + np.ndarray: ``(a + b) / 2``. + + Source: + Stillstanding study, pose_motion.py (shoulder-midpoint torso + micromotion) (Jensenius). + """ + return (np.asarray(a, dtype=np.float64) + np.asarray(b, dtype=np.float64)) / 2.0 + + +def limb_speed_from_landmarks( + xy: np.ndarray, + confidence: np.ndarray | None, + fps: float, + conf_gate: float = 0.5, + merge: str | None = "max_lr", + smooth_taps: int = 3) -> np.ndarray: + """ + Confidence-gated image-plane speed of one or more candidate limbs. + + For each candidate limb (e.g. the left and right wrist), frames whose + landmark confidence/visibility falls below ``conf_gate`` are masked out + (NaN), and the limb speed is formed as the central-difference magnitude of + the pixel path (px/s). Candidate limbs are then merged by element-wise + maximum — so that motion of *either* limb registers, mirroring the + bilateral merge used for inertial hand data — and lightly smoothed with a + short NaN-aware moving average. Peaks of the resulting signal mark, e.g., + strike downstrokes of the striking wrist. + + Caveats (from the cymbal study): these are 2D apparent kinematics from a + single camera — motion toward/away from the lens is foreshortened and pixel + speed is not metric speed. Moreover, a limb-speed peak marks *maximum + downstroke speed*, which systematically precedes the contact/arrest that an + audio onset or an acceleration peak registers; account for this bias when + comparing event times across modalities. + + Peak-picking on the returned signal is left to the caller (a general + adaptive peak-picker, ``pick_peaks``, is provided by the sibling + core-signal-methods PR in ``musicalgestures._peaks``; the cymbal study used + a relative threshold of 0.4 x the take's peak with a 0.2 s minimum + interval). + + Args: + xy (np.ndarray): Pixel positions, shape (F, L, 2) for L candidate limbs + or (F, 2) for a single limb. + confidence (np.ndarray, optional): Per-frame landmark confidence + (MediaPipe ``visibility``), shape (F, L) or (F,). Pass None to skip + confidence gating. + fps (float): Frame rate of the trajectory (Hz). + conf_gate (float, optional): Frames with confidence below this value + are masked (NaN) before differentiation. Defaults to 0.5. + merge (str, optional): ``"max_lr"`` (or ``"max"``) merges the candidate + limbs by element-wise (NaN-aware) maximum; None returns per-limb + speeds. Defaults to "max_lr". + smooth_taps (int, optional): Length of the NaN-aware moving-average + smoother applied after merging. Use 0 or 1 to disable. Defaults + to 3. + + Returns: + np.ndarray: Speed in px/s: shape (F,) when merged, else (F, L). NaN + where the position (or a central-difference neighbour) is masked + or missing. + + Source: + Cymbal-comparison study, markerless striking-wrist speed + (reimplemented from the paper's method description; defaults are the + paper's provisional values) (Jensenius). + """ + xy = np.asarray(xy, dtype=np.float64) + single = xy.ndim == 2 + if single: + xy = xy[:, None, :] + if xy.ndim != 3 or xy.shape[-1] != 2: + raise ValueError(f"xy must have shape (F, L, 2) or (F, 2), got {xy.shape}") + xy = xy.copy() + + if confidence is not None: + conf = np.asarray(confidence, dtype=np.float64) + if conf.ndim == 1: + conf = conf[:, None] + if conf.shape != xy.shape[:2]: + raise ValueError( + f"confidence shape {conf.shape} does not match xy frames/limbs {xy.shape[:2]}") + with np.errstate(invalid="ignore"): + xy[~(conf >= conf_gate)] = np.nan + + if len(xy) < 2: + speed = np.full(xy.shape[:2], np.nan) + else: + # Central differences (one-sided at the ends); NaN wherever a needed + # neighbour is masked/missing. + d = np.gradient(xy, axis=0) + speed = np.hypot(d[..., 0], d[..., 1]) * float(fps) + + if merge is not None: + if merge not in ("max_lr", "max"): + raise ValueError(f"merge must be 'max_lr', 'max' or None, got {merge!r}") + with warnings.catch_warnings(): + warnings.simplefilter("ignore", category=RuntimeWarning) # all-NaN frames + speed = np.nanmax(speed, axis=1) + elif single: + speed = speed[:, 0] + + if smooth_taps and smooth_taps > 1: + speed = _nan_moving_average(speed, int(smooth_taps)) + return speed + + +def impact_events( + pos_by_point: np.ndarray, + fps: float, + rel_thresh: float = 0.12, + min_interval_s: float = 0.10) -> dict: + """ + Detect candidate impact events from point trajectories via acceleration peaks. + + Each candidate point's position (2D or 3D; e.g. the two hand points of a + mocap model, or the two wrist landmarks) is differentiated twice with + central differences to obtain its acceleration vector, the vector magnitude + is taken, and the points are merged by element-wise maximum so that a + strike by *either* hand registers (bilateral max). Impacts are then + peak-picked on the merged acceleration magnitude with a relative threshold + of ``rel_thresh`` x the signal's maximum and a minimum inter-impact + interval of ``min_interval_s``. + + The threshold parameters are taken directly (the small relative-threshold + peak picker is implemented inline here); a general adaptive peak-picker, + ``pick_peaks``, is provided by the sibling core-signal-methods PR in + ``musicalgestures._peaks``. The defaults (0.12 x peak, 100 ms) are the + cymbal study's provisional values for 120 Hz mocap hand data and should be + tuned per dataset. Note the study's caveat: double-differentiating + (model-reconstructed) positions is noisy and also responds to the + backswing, not only the collision — treat the detected peaks as *candidate* + impacts and validate against another modality (e.g. audio onsets) where + possible. For whole-image visual impact detection from video (no + landmarks), see ``MgVideo.impacts()`` instead. + + Args: + pos_by_point (np.ndarray): Point positions, shape (F, P, D) for P + candidate points in D spatial dimensions (2 or 3), or (F, D) for a + single point. Units are the caller's (m or px); NaNs (dropouts) + propagate into the acceleration and are never picked as peaks. + fps (float): Sampling rate of the trajectories (Hz). + rel_thresh (float, optional): Relative peak threshold as a fraction of + the merged acceleration magnitude's maximum. Defaults to 0.12. + min_interval_s (float, optional): Minimum interval between detected + impacts in seconds (stronger peaks win). Defaults to 0.10. + + Returns: + dict: A dictionary with keys: + + - ``index`` (np.ndarray of int): Sample indices of the detected + impacts, ascending. + - ``time`` (np.ndarray): Impact times in seconds (``index / fps``). + - ``magnitude`` (np.ndarray): Merged acceleration magnitude at each + impact (position-units/s^2). + - ``accel`` (np.ndarray, shape (F,)): The full merged acceleration- + magnitude signal (for plotting/inspection). + + Source: + Cymbal-comparison study, kinematic impact detection from Xsens hand + points (reimplemented from the paper's method description; defaults + are the paper's provisional values) (Jensenius). + """ + pos = np.asarray(pos_by_point, dtype=np.float64) + if pos.ndim == 2: + pos = pos[:, None, :] + if pos.ndim != 3: + raise ValueError( + f"pos_by_point must have shape (F, P, D) or (F, D), got {pos.shape}") + + n = len(pos) + if n < 3: + accel = np.full(n, np.nan) + else: + vel = np.gradient(pos, axis=0) * float(fps) + acc = np.gradient(vel, axis=0) * float(fps) + mag = np.linalg.norm(acc, axis=2) # (F, P) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", category=RuntimeWarning) # all-NaN frames + accel = np.nanmax(mag, axis=1) # bilateral max across points + + min_dist = max(1, int(round(min_interval_s * float(fps)))) + idx = _pick_relative_peaks(accel, rel_thresh, min_dist) + return { + "index": idx, + "time": idx / float(fps), + "magnitude": accel[idx] if len(idx) else np.array([], dtype=np.float64), + "accel": accel, + } + + +def _nan_moving_average(x: np.ndarray, taps: int) -> np.ndarray: + """NaN-aware centred moving average along the first axis. + + Averages the finite values inside each window; a sample is NaN only when + its whole window is NaN. ``taps`` should be odd (an even value is widened + by one to stay centred). + """ + x = np.asarray(x, dtype=np.float64) + if taps % 2 == 0: + taps += 1 + half = taps // 2 + stack = np.full((taps,) + x.shape, np.nan) + for k in range(-half, half + 1): + src = slice(max(0, -k), x.shape[0] - max(0, k)) + dst = slice(max(0, k), x.shape[0] - max(0, -k)) + stack[k + half][dst] = x[src] + with warnings.catch_warnings(): + warnings.simplefilter("ignore", category=RuntimeWarning) # all-NaN windows + return np.nanmean(stack, axis=0) + + +def _pick_relative_peaks(signal: np.ndarray, rel_thresh: float, min_dist: int) -> np.ndarray: + """Pick local maxima above ``rel_thresh`` x max, at least ``min_dist`` samples apart. + + Minimal inline relative-threshold picker (stronger peaks win ties within + ``min_dist``); NaN samples are never candidates. Kept private — the + general, feature-complete picker lives in ``musicalgestures._peaks`` + (sibling PR). + """ + s = np.asarray(signal, dtype=np.float64) + if len(s) < 3 or not np.isfinite(s).any(): + return np.array([], dtype=np.intp) + peak = np.nanmax(s) + if not peak > 0: + return np.array([], dtype=np.intp) + thr = rel_thresh * peak + interior = np.arange(1, len(s) - 1) + mid, prev, nxt = s[1:-1], s[:-2], s[2:] + cand = interior[(mid >= thr) & (mid > prev) & (mid >= nxt)] + if len(cand) == 0: + return cand.astype(np.intp) + kept: list[int] = [] + for i in cand[np.argsort(s[cand], kind="stable")[::-1]]: + if all(abs(int(i) - j) >= min_dist for j in kept): + kept.append(int(i)) + return np.array(sorted(kept), dtype=np.intp) From 377e3c9beb57af3e3e28489c85bf8fd43eddf835 Mon Sep 17 00:00:00 2001 From: Alexander Refsum Jensenius Date: Thu, 16 Jul 2026 09:43:39 +0200 Subject: [PATCH 2/3] Add tests for _posetools Synthetic ground-truth tests (always run) for midpoint, limb_speed_from_landmarks and impact_events: known constant-velocity speeds, confidence gating, bilateral max merges, Gaussian-bump acceleration impacts recovered after double integration, relative threshold and minimum-interval behaviour, NaN dropout handling. Import-safety is verified in a subprocess with mediapipe blocked (package import and numpy-only helpers work; extract_pose_landmarks raises an ImportError naming the [pose] extra). The extractor itself is covered by an integration test on the bundled dancer example video that skips when mediapipe or the model file is unavailable. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_0147io2kMk8M6jcNFNt1G96m --- tests/test_posetools.py | 331 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 331 insertions(+) create mode 100644 tests/test_posetools.py diff --git a/tests/test_posetools.py b/tests/test_posetools.py new file mode 100644 index 0000000..f7f5519 --- /dev/null +++ b/tests/test_posetools.py @@ -0,0 +1,331 @@ +"""Tests for musicalgestures._posetools (landmark-trajectory pose tools). + +The numpy-only helpers (midpoint, limb_speed_from_landmarks, impact_events) +are tested against synthetic trajectories with known ground truth and always +run. extract_pose_landmarks needs MediaPipe (optional dependency) plus a real +human video, so it is an integration test that skips when either is missing. +""" + +import os +import subprocess +import sys + +import numpy as np +import pytest + +import musicalgestures +from musicalgestures._posetools import ( + extract_pose_landmarks, + midpoint, + limb_speed_from_landmarks, + impact_events, + _pick_relative_peaks, +) + + +def _gaussian_bump(n, center, sigma, amplitude): + t = np.arange(n, dtype=float) + return amplitude * np.exp(-0.5 * ((t - center) / sigma) ** 2) + + +def _position_from_accel(accel, fps): + """Double-integrate a 1-D acceleration signal into a position signal.""" + vel = np.cumsum(accel) / fps + return np.cumsum(vel) / fps + + +# --------------------------------------------------------------------------- +# midpoint +# --------------------------------------------------------------------------- + +class TestMidpoint: + def test_basic(self): + a = np.array([[0.0, 0.0], [2.0, 2.0]]) + b = np.array([[2.0, 4.0], [4.0, 6.0]]) + np.testing.assert_allclose(midpoint(a, b), [[1.0, 2.0], [3.0, 4.0]]) + + def test_nan_propagates(self): + a = np.array([[0.0, np.nan], [2.0, 2.0]]) + b = np.array([[2.0, 4.0], [4.0, 6.0]]) + m = midpoint(a, b) + assert np.isnan(m[0, 1]) + np.testing.assert_allclose(m[1], [3.0, 4.0]) + + def test_shoulder_midpoint_from_landmark_array(self): + # (F, 33, 3) landmark array as returned by extract_pose_landmarks + lm = np.zeros((4, 33, 3)) + lm[:, 11, :2] = [100.0, 50.0] # left shoulder + lm[:, 12, :2] = [200.0, 70.0] # right shoulder + m = midpoint(lm[:, 11, :2], lm[:, 12, :2]) + np.testing.assert_allclose(m, np.tile([150.0, 60.0], (4, 1))) + + +# --------------------------------------------------------------------------- +# limb_speed_from_landmarks +# --------------------------------------------------------------------------- + +class TestLimbSpeed: + def test_constant_velocity_speed(self): + # One limb moving at (3, 4) px/frame at 50 fps -> speed 5 px/frame = 250 px/s + fps = 50.0 + t = np.arange(100, dtype=float) + xy = np.stack([3.0 * t, 4.0 * t], axis=1) + conf = np.ones(100) + speed = limb_speed_from_landmarks(xy, conf, fps) + assert speed.shape == (100,) + np.testing.assert_allclose(speed, 250.0) + + def test_confidence_gate_masks_frames(self): + fps = 50.0 + t = np.arange(100, dtype=float) + xy = np.stack([3.0 * t, 4.0 * t], axis=1) + conf = np.ones(100) + conf[40:51] = 0.2 # below the default 0.5 gate + speed = limb_speed_from_landmarks(xy, conf, fps) + assert np.isnan(speed[45]) + # away from the gated region (and its central-difference/smoothing + # neighbourhood) the speed is intact + np.testing.assert_allclose(speed[5:35], 250.0) + np.testing.assert_allclose(speed[60:95], 250.0) + + def test_bilateral_max_merge(self): + # Left limb moves at 5 px/frame in the first half, right limb at + # 10 px/frame in the second half; merged speed tracks the faster limb. + fps = 30.0 + n = 120 + left = np.zeros((n, 2)) + right = np.zeros((n, 2)) + left[:60, 0] = 5.0 * np.arange(60) + left[60:, 0] = left[59, 0] + right[:60] = 0.0 + right[60:, 1] = 10.0 * np.arange(60) + xy = np.stack([left, right], axis=1) # (F, 2, 2) + conf = np.ones((n, 2)) + merged = limb_speed_from_landmarks(xy, conf, fps, smooth_taps=0) + per_limb = limb_speed_from_landmarks(xy, conf, fps, merge=None, smooth_taps=0) + assert per_limb.shape == (n, 2) + # interior of each half: merged = max of the two limbs + np.testing.assert_allclose(merged[10:50], 5.0 * fps) + np.testing.assert_allclose(merged[70:110], 10.0 * fps) + np.testing.assert_allclose(merged, np.nanmax(per_limb, axis=1)) + + def test_speed_peak_location(self): + # Velocity is a Gaussian bump centred at frame 100: the speed peak + # must sit at the bump centre. + fps = 60.0 + n = 200 + vel = _gaussian_bump(n, center=100, sigma=6.0, amplitude=8.0) # px/frame + xy = np.stack([np.cumsum(vel), np.zeros(n)], axis=1) + speed = limb_speed_from_landmarks(xy, np.ones(n), fps) + assert abs(int(np.nanargmax(speed)) - 100) <= 1 + # peak speed close to the designed maximum (px/frame * fps) + assert speed[np.nanargmax(speed)] == pytest.approx(8.0 * fps, rel=0.05) + + def test_no_confidence_given(self): + fps = 25.0 + t = np.arange(50, dtype=float) + xy = np.stack([2.0 * t, 0.0 * t], axis=1) + speed = limb_speed_from_landmarks(xy, None, fps) + np.testing.assert_allclose(speed, 50.0) + + def test_all_limbs_gated_gives_nan(self): + xy = np.zeros((20, 2, 2)) + conf = np.zeros((20, 2)) + speed = limb_speed_from_landmarks(xy, conf, 30.0) + assert speed.shape == (20,) + assert np.isnan(speed).all() + + def test_bad_shapes_raise(self): + with pytest.raises(ValueError): + limb_speed_from_landmarks(np.zeros((10, 3)), None, 30.0) + with pytest.raises(ValueError): + limb_speed_from_landmarks(np.zeros((10, 2, 2)), np.ones((9, 2)), 30.0) + with pytest.raises(ValueError): + limb_speed_from_landmarks(np.zeros((10, 2)), None, 30.0, merge="mean") + + +# --------------------------------------------------------------------------- +# impact_events +# --------------------------------------------------------------------------- + +class TestImpactEvents: + fps = 120.0 + + def _pos_from_bumps(self, n, bumps): + """1-point 2-D position whose x-acceleration has Gaussian bumps. + + bumps: list of (center_index, amplitude). + """ + accel = np.zeros(n) + for c, a in bumps: + accel += _gaussian_bump(n, c, sigma=2.0, amplitude=a) + x = _position_from_accel(accel, self.fps) + return np.stack([x, np.zeros(n)], axis=1) # (F, 2) + + def test_detects_known_impacts(self): + pos = self._pos_from_bumps(600, [(100, 10.0), (250, 6.0), (400, 4.0)]) + ev = impact_events(pos, self.fps) + assert len(ev["index"]) == 3 + for found, expected in zip(ev["index"], (100, 250, 400)): + assert abs(found - expected) <= 2 + np.testing.assert_allclose(ev["time"], ev["index"] / self.fps) + assert len(ev["accel"]) == 600 + # strongest impact carries the largest magnitude + assert np.argmax(ev["magnitude"]) == 0 + + def test_bilateral_max_across_points(self): + # A strike by either hand registers: one bump per point. + left = self._pos_from_bumps(600, [(150, 8.0)]) + right = self._pos_from_bumps(600, [(350, 8.0)]) + pos = np.stack([left, right], axis=1) # (F, 2, 2) + ev = impact_events(pos, self.fps) + assert len(ev["index"]) == 2 + assert abs(ev["index"][0] - 150) <= 2 + assert abs(ev["index"][1] - 350) <= 2 + + def test_relative_threshold_excludes_weak_peaks(self): + pos = self._pos_from_bumps(600, [(100, 10.0), (300, 0.5)]) # 5% of max + ev = impact_events(pos, self.fps, rel_thresh=0.12) + assert len(ev["index"]) == 1 + assert abs(ev["index"][0] - 100) <= 2 + # a permissive threshold keeps the weak one too + ev2 = impact_events(pos, self.fps, rel_thresh=0.02) + assert len(ev2["index"]) == 2 + + def test_min_interval_keeps_stronger_peak(self): + # Two impacts 6 samples (50 ms) apart; min_interval 100 ms keeps only + # the stronger. + pos = self._pos_from_bumps(600, [(200, 10.0), (206, 6.0)]) + ev = impact_events(pos, self.fps, min_interval_s=0.10) + close = [i for i in ev["index"] if 190 <= i <= 216] + assert len(close) == 1 + assert abs(close[0] - 200) <= 2 + # with a permissive interval both survive + ev2 = impact_events(pos, self.fps, min_interval_s=0.02) + close2 = [i for i in ev2["index"] if 190 <= i <= 216] + assert len(close2) == 2 + + def test_3d_positions(self): + n = 600 + accel = _gaussian_bump(n, 300, sigma=2.0, amplitude=9.0) + x = _position_from_accel(accel, self.fps) + pos = np.stack([x, x, x], axis=1)[:, None, :] # (F, 1, 3) + ev = impact_events(pos, self.fps) + assert len(ev["index"]) == 1 + assert abs(ev["index"][0] - 300) <= 2 + + def test_nan_dropout_region(self): + pos = self._pos_from_bumps(600, [(100, 10.0), (400, 8.0)]) + pos[200:250] = np.nan + ev = impact_events(pos, self.fps) + found = sorted(ev["index"]) + assert any(abs(i - 100) <= 2 for i in found) + assert any(abs(i - 400) <= 2 for i in found) + assert not any(200 <= i < 250 for i in found) + + def test_empty_and_flat_signals(self): + ev = impact_events(np.zeros((2, 2)), self.fps) + assert len(ev["index"]) == 0 + ev = impact_events(np.ones((100, 2)), self.fps) # no motion at all + assert len(ev["index"]) == 0 + + def test_bad_shape_raises(self): + with pytest.raises(ValueError): + impact_events(np.zeros((10, 2, 2, 2)), self.fps) + + +class TestPickRelativePeaks: + def test_plateau_and_order(self): + s = np.array([0.0, 1.0, 1.0, 0.0, 5.0, 0.0, 2.0, 0.0]) + idx = _pick_relative_peaks(s, rel_thresh=0.1, min_dist=1) + assert list(idx) == [1, 4, 6] # plateau picked once, ascending order + + def test_nan_never_a_peak(self): + s = np.array([0.0, np.nan, 0.0, 3.0, 0.0]) + idx = _pick_relative_peaks(s, rel_thresh=0.1, min_dist=1) + assert list(idx) == [3] + + +# --------------------------------------------------------------------------- +# extract_pose_landmarks +# --------------------------------------------------------------------------- + +def test_import_safe_without_mediapipe(): + """The package (and the numpy-only helpers) must work without mediapipe, + and extract_pose_landmarks must raise a clear ImportError naming the + [pose] extra.""" + code = "\n".join([ + "import sys", + "sys.modules['mediapipe'] = None # simulate mediapipe not installed", + "import numpy as np", + "import musicalgestures", + "m = musicalgestures.midpoint(np.array([0.0, 0.0]), np.array([2.0, 4.0]))", + "assert m.tolist() == [1.0, 2.0]", + "try:", + " musicalgestures.extract_pose_landmarks('does_not_matter.mp4')", + "except ImportError as exc:", + " assert 'musicalgestures[pose]' in str(exc), str(exc)", + "else:", + " raise AssertionError('expected ImportError without mediapipe')", + "print('IMPORT_SAFE_OK')", + ]) + out = subprocess.run([sys.executable, "-c", code], capture_output=True, text=True) + assert out.returncode == 0, out.stderr + assert "IMPORT_SAFE_OK" in out.stdout + + +def _mediapipe_available(): + try: + import mediapipe # noqa: F401 + return True + except Exception: + return False + + +EXAMPLE_VIDEO = os.path.join( + os.path.dirname(musicalgestures.__file__), "examples", "dancer.avi") + + +@pytest.mark.skipif(not _mediapipe_available(), reason="mediapipe not installed") +@pytest.mark.skipif(not os.path.exists(EXAMPLE_VIDEO), reason="example video missing") +def test_extract_pose_landmarks_integration(tmp_path): + """Integration test on the bundled dancer example video (real human + motion). Skips when the MediaPipe model file cannot be obtained.""" + from musicalgestures._exceptions import MgDependencyError + + csv_path = str(tmp_path / "dancer_landmarks.csv") + try: + res = extract_pose_landmarks( + EXAMPLE_VIDEO, fps=5, width=320, world_landmarks=True, + max_frames=15, target_name=csv_path, verbose=False) + except MgDependencyError as exc: # model download failed (offline) + pytest.skip(f"MediaPipe model unavailable: {exc}") + + n = len(res["time"]) + assert 0 < n <= 15 + assert res["landmarks"].shape == (n, 33, 3) + assert res["world"].shape == (n, 33, 3) + assert res["detected"].shape == (n,) + assert res["width"] == 320 + assert res["fps"] == 5 + assert len(res["names"]) == 33 + # the dancer is plainly visible: expect a solid detection rate + assert res["detection_rate"] > 0.5 + # detected frames have finite pixel coordinates within the analysis frame, + # undetected frames are all-NaN + det = res["detected"] + assert np.isfinite(res["landmarks"][det][:, :, :2]).all() + if (~det).any(): + assert np.isnan(res["landmarks"][~det]).all() + x = res["landmarks"][det][:, :, 0] + assert (x > -res["width"]).all() and (x < 2 * res["width"]).all() + # timestamps follow the analysis rate + np.testing.assert_allclose(np.diff(res["time"]), 1.0 / 5.0) + # tidy CSV written with one column per landmark coordinate + assert os.path.exists(csv_path) + with open(csv_path) as fh: + header = fh.readline().strip().split(",") + assert header[0] == "time" + assert "left_wrist_x" in header and "left_wrist_v" in header + assert "left_wrist_wz" in header # world landmarks requested + assert len(header) == 1 + 33 * 3 + 33 * 3 From eabd9e74fe31b0c774b8724db4f491e829601238 Mon Sep 17 00:00:00 2001 From: Alexander Refsum Jensenius Date: Thu, 16 Jul 2026 09:56:21 +0200 Subject: [PATCH 3/3] fix: PR3 review findings (legacy-path markers, stderr drain, validation consistency) Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_0147io2kMk8M6jcNFNt1G96m --- musicalgestures/_posetools.py | 49 ++++++++++++++++++++++++++++++++--- tests/test_posetools.py | 8 ++++++ 2 files changed, 54 insertions(+), 3 deletions(-) diff --git a/musicalgestures/_posetools.py b/musicalgestures/_posetools.py index 108d711..9303cfd 100644 --- a/musicalgestures/_posetools.py +++ b/musicalgestures/_posetools.py @@ -31,6 +31,7 @@ import os import subprocess +import threading import warnings import numpy as np @@ -130,6 +131,10 @@ def extract_pose_landmarks( (mp_extract_westney.py, pose_motion.py) and Westney-comparisons (concert_mediapipe.py, reh_pose.py, a1_labstage.py) (Jensenius). """ + if model_complexity not in (0, 1, 2): + raise ValueError( + f"model_complexity must be 0, 1 or 2, got {model_complexity!r}") + try: import mediapipe as mp except ImportError as exc: @@ -168,6 +173,22 @@ def extract_pose_landmarks( # The study scripts were written against the legacy Solutions API # (mediapipe 0.10.14); newer 0.10.x wheels (e.g. 0.10.35, as used in the # cymbal study) removed it in favour of the Tasks API. Support both. + # + # `use_solutions` prefers the Solutions API whenever it is present: this + # is purely for fidelity to the original study pipeline (the papers' + # numbers were produced with Solutions, not Tasks), not because Solutions + # is otherwise preferable. + # + # IMPORTANT: the Solutions branch below (`if use_solutions:`) is a + # faithful port of the study scripts' mediapipe<=0.10.14 Solutions code, + # but it is UNTESTABLE on modern wheels — mediapipe>=0.10.26 removed + # `mp.solutions` entirely, so no environment this project can currently + # install exercises this branch (it is only reachable with an old, pinned + # mediapipe wheel). Keep it faithful to the source scripts rather than + # "improving" it blind. + # TODO: once the legacy Solutions-API mediapipe family (<=0.10.14) is + # fully retired (no supported environment can install it any more), drop + # this preference and the Solutions branch, and always use the Tasks API. use_solutions = hasattr(mp, "solutions") def _read_result(lms, wlms): @@ -184,6 +205,20 @@ def _read_result(lms, wlms): frame_bytes = w * h * 3 stopped_early = False proc = subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE) + # Drain FFmpeg's stderr on a background thread as it is produced, rather + # than only reading it in the `finally` block below: the stdout-reading + # loop can run far longer than the OS pipe buffer takes to fill (a few + # dozen KB), and if FFmpeg blocks on a full stderr pipe while we are only + # consuming stdout, decoding stalls. The drained chunks are joined at the + # end to preserve the existing error-reporting behaviour. + stderr_chunks: list[bytes] = [] + + def _drain_stderr(): + for chunk in iter(lambda: proc.stderr.read(4096), b""): + stderr_chunks.append(chunk) + + stderr_thread = threading.Thread(target=_drain_stderr, daemon=True) + stderr_thread.start() try: with _suppress_native_stderr(quiet): if use_solutions: @@ -206,6 +241,8 @@ def _process(frame_rgb, t): # Tasks API: reuse the model download/cache logic of the # per-frame estimator so both pose workflows share one model # file in musicalgestures/models/. + # TODO: promote _get_model_path to a public helper in + # _pose_estimator; this is a private cross-module call. from musicalgestures._pose_estimator import MediaPipePoseEstimator model_path = MediaPipePoseEstimator( model_complexity=model_complexity)._get_model_path() @@ -257,9 +294,10 @@ def _process(frame_rgb, t): close_backend() finally: proc.stdout.close() - err = proc.stderr.read().decode(errors="replace").strip() - proc.stderr.close() proc.wait() + stderr_thread.join() + proc.stderr.close() + err = b"".join(stderr_chunks).decode(errors="replace").strip() # Stopping at max_frames breaks FFmpeg's output pipe on purpose, so # suppress its resulting broken-pipe complaints. if err and not stopped_early: @@ -379,7 +417,12 @@ def limb_speed_from_landmarks( speeds. Defaults to "max_lr". smooth_taps (int, optional): Length of the NaN-aware moving-average smoother applied after merging. Use 0 or 1 to disable. Defaults - to 3. + to 3. With ``smooth_taps > 1``, samples right at the edge of a + confidence-gated (NaN) region can be partially reconstructed: + the NaN-aware average only requires *some* finite values inside + its window, so an edge sample whose window straddles both valid + and gated frames is averaged from the valid ones rather than + staying NaN. Returns: np.ndarray: Speed in px/s: shape (F,) when merged, else (F, L). NaN diff --git a/tests/test_posetools.py b/tests/test_posetools.py index f7f5519..6b8fe2d 100644 --- a/tests/test_posetools.py +++ b/tests/test_posetools.py @@ -274,6 +274,14 @@ def test_import_safe_without_mediapipe(): assert "IMPORT_SAFE_OK" in out.stdout +def test_invalid_model_complexity_raises(): + """model_complexity must be 0, 1 or 2; validated up front so both the + Solutions and Tasks backends behave identically (and mediapipe need not + be installed for this to raise).""" + with pytest.raises(ValueError): + extract_pose_landmarks("does_not_matter.mp4", model_complexity=3) + + def _mediapipe_available(): try: import mediapipe # noqa: F401