From 62543113680507570a5ef79a81f30ae4b7ef124b Mon Sep 17 00:00:00 2001 From: Alexander Refsum Jensenius Date: Wed, 29 Jul 2026 09:27:54 +0200 Subject: [PATCH] Depend on micromotion instead of carrying our own quantity-of-motion code MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit _qom, _mocap, _posture and _physio move to the micromotion package and become re-export shims. Behaviour is unchanged: every one of the 53 tests in test_qom/test_mocap/test_posture/test_physio passes against micromotion unmodified, and band_limited_qom still returns 5.675 mm/s on the reference 200 Hz optical recording. Why this direction. These functions are credited in their own docstrings to the Still Standing study, which was independently extracting the same measure into micromotion — so the project's central measure was about to exist in two disagreeing copies inside one lab. micromotion needs only numpy, scipy and pandas, so depending on it costs nothing, whereas the reverse would make anyone analysing accelerometer data install a computer-vision stack. It matches the arrow MGT already has to ambiscape. Two duplicate implementations were collapsed onto one code path in micromotion: dfa (float wrapper over the dict-returning version; they agreed to 1.2 per cent on Brownian motion) and sample_entropy (agreed to 0.2 per cent). One item is deliberately left alone: axial_rayleigh. Ours returns R = 0.9985 where micromotion's circular.rayleigh_axial returns 0.9222 for the same bimodal input, and only the latter matches the textbook resultant length of the doubled angles. Ours is kept as-is here because it is the one under test; the discrepancy needs resolving on its own rather than inside this move. --- musicalgestures/_mocap.py | 230 ++----------- musicalgestures/_physio.py | 153 +-------- musicalgestures/_posture.py | 667 +++--------------------------------- musicalgestures/_qom.py | 412 ++-------------------- pyproject.toml | 1 + 5 files changed, 112 insertions(+), 1351 deletions(-) diff --git a/musicalgestures/_mocap.py b/musicalgestures/_mocap.py index 5085fe3..dcd8b36 100644 --- a/musicalgestures/_mocap.py +++ b/musicalgestures/_mocap.py @@ -1,214 +1,26 @@ -""" -Motion-capture I/O and cross-modality utilities. - -Pure numpy/scipy helpers ported from the "still standing" and -Westney-comparisons studies: +"""Re-exported from the ``micromotion`` package. -* :func:`read_qtm_tsv` -- a single robust reader for Qualisys Track Manager - (QTM) TSV exports, consolidating the four/five near-duplicate loaders that - were copy-pasted across the study scripts. -* :func:`compare_modality_envelopes` -- resample two motion envelopes onto a - common per-second grid and correlate them (e.g. video-pose vs mocap - validation). -* :func:`dominant_frequency` -- the dominant spectral peak of a signal within - a band, via Welch. +These functions used to live here. They were moved to ``micromotion`` on 2026-07-29 so that +one implementation of quantity of motion exists rather than two, and MGT now depends on that +package instead of carrying its own copy. Behaviour is unchanged: this module's tests pass +against ``micromotion`` unmodified. -.. note:: - :func:`compare_modality_envelopes` deliberately takes *precomputed* 1-D - motion envelopes rather than computing quantity-of-motion internally, so - this module stays independent of the QoM machinery. The natural producer - of such envelopes is ``musicalgestures._qom.band_limited_qom`` followed by - a per-second binning (``envelope`` / ``bin_series``), arriving in a - sibling PR. +The dependency points this way round on purpose. ``micromotion`` needs only numpy, scipy and +pandas, so someone analysing accelerometer data does not have to install a computer-vision +stack; MGT already depends on ``ambiscape`` the same way, and neither of those packages +imports MGT. -Source: still standing study and Westney-comparisons study (Jensenius). +Import from ``micromotion`` directly in new code. """ -import numpy as np - - -def read_qtm_tsv(path): - """ - Read a Qualisys Track Manager (QTM) TSV motion-capture export. - - Consolidates the several near-duplicate loaders used across the studies - into one robust reader. It locates the ``MARKER_NAMES`` header row to - recover marker labels, autodetects where the numeric data block starts - (the first row whose first field parses as a float), drops a trailing - all-empty column produced by a trailing tab, converts exact-zero XYZ - triples (Qualisys gap fills) to ``NaN``, and falls back from UTF-8 to - latin-1 encoding. When a ``FREQUENCY`` header field is present the frame - rate is returned as well. - - Source: still standing study and Westney-comparisons study (Jensenius) -- - unifies the ``load_qtm`` variants in the balance/dynamics/circular/ - spatial-range reports and the latin-1 variant in ``compare_mp_mocap``. - - Args: - path (str): Path to the ``.tsv`` file. - - Returns: - tuple: ``(marker_names, data, fs)`` where ``marker_names`` is a list - of ``M`` strings (empty if no header was found), ``data`` is a - float array of shape ``(T, M, 3)`` with gaps as ``NaN``, and - ``fs`` is the frame rate in Hz or ``None`` if not derivable from - the header. - """ - import io - - def _read_lines(enc): - with io.open(path, encoding=enc) as fh: - return fh.readlines() - - try: - lines = _read_lines("utf-8") - except (UnicodeDecodeError, UnicodeError): - lines = _read_lines("latin-1") - - marker_names = [] - fs = None - data_start = None - for i, ln in enumerate(lines): - parts = ln.rstrip("\n").split("\t") - key = parts[0].strip().upper() - if key == "MARKER_NAMES": - marker_names = [p.strip() for p in parts[1:] if p.strip()] - continue - if key == "FREQUENCY" and len(parts) > 1: - try: - fs = float(parts[1]) - except ValueError: - pass - continue - # first row whose leading field is numeric marks the data block - try: - float(parts[0]) - data_start = i - break - except ValueError: - continue - - if data_start is None: - raise ValueError(f"No numeric data block found in {path!r}") - - rows = [] - for ln in lines[data_start:]: - parts = ln.rstrip("\n").split("\t") - try: - rows.append([float(p) for p in parts if p != ""]) - except ValueError: - continue - if not rows: - raise ValueError(f"No parseable data rows in {path!r}") - - ncol = min(len(r) for r in rows) - arr = np.array([r[:ncol] for r in rows], dtype=float) - - # A QTM data row may carry a leading time/frame column(s); the marker - # block is the trailing 3*M columns. Prefer the marker count from the - # header; otherwise infer the largest multiple of 3. - if marker_names: - M = len(marker_names) - else: - M = ncol // 3 - offset = arr.shape[1] - 3 * M - if offset < 0: - # header lists more markers than columns present; fall back - M = arr.shape[1] // 3 - offset = arr.shape[1] - 3 * M - marker_names = marker_names[:M] - block = arr[:, offset:offset + 3 * M] - data = block.reshape(block.shape[0], M, 3) - - # exact-zero triples are Qualisys gap fills -> NaN - gap = np.all(data == 0, axis=2) - data[gap] = np.nan - - return marker_names, data, fs - - -def dominant_frequency(x, fs, band=(0.3, 4.0)): - """ - Dominant frequency of a signal within a band, via a Welch spectrum. - - Returns the frequency of the largest Welch power-spectral-density peak - inside ``band`` -- e.g. the dominant oscillation rate of a body-part - speed or vertical-position signal. - - Source: Westney-comparisons study (Jensenius), extended motion-feature - analysis (motion dominant frequency of vertical trunk position). - - Args: - x (np.ndarray): 1-D input signal. - fs (float): Sampling rate in Hz. - band (tuple, optional): ``(low, high)`` search band in Hz. Defaults - to ``(0.3, 4.0)``. - - Returns: - float: The dominant frequency in Hz, or ``nan`` if the band is empty. - """ - from scipy.signal import welch - - x = np.asarray(x, dtype=float) - x = x[np.isfinite(x)] - if len(x) < 8: - return np.nan - lo, hi = band - f, P = welch(x - x.mean(), fs, nperseg=min(2048, len(x))) - mask = (f >= lo) & (f <= hi) - if not mask.any(): - return np.nan - return float(f[mask][np.argmax(P[mask])]) - - -def compare_modality_envelopes(env_a, env_b, fs_a, fs_b): - """ - Correlate two motion envelopes after resampling to a common grid. - - Resamples both 1-D envelopes onto a shared one-sample-per-second grid - (by averaging within each second), truncates to the common length, and - returns the Pearson correlation -- the video-vs-mocap (or view-vs-view) - agreement measure. Both inputs are treated as already-computed motion - envelopes (e.g. per-frame band-limited quantity-of-motion), keeping this - function decoupled from the QoM computation itself. The per-second binning - uses an integer-rounded step, so non-integer frame rates (e.g. 29.97 fps) - drift slightly over long signals; this function is intended for validation - rather than precise alignment. - - Source: still standing / Westney-comparisons study (Jensenius), - MediaPipe-vs-mocap validation (``compare_mp_mocap``). - - Args: - env_a (np.ndarray): First 1-D motion envelope. - env_b (np.ndarray): Second 1-D motion envelope. - fs_a (float): Sampling rate of ``env_a`` in Hz. - fs_b (float): Sampling rate of ``env_b`` in Hz. - - Returns: - dict: ``{"r", "n"}`` where ``r`` is the Pearson correlation of the - two per-second envelopes and ``n`` the number of common seconds - (``r`` is ``nan`` if fewer than three overlapping seconds or if - either resampled envelope is constant). - """ - def _per_second(env, fs): - env = np.asarray(env, dtype=float) - nsec = int(len(env) / fs) - if nsec < 1: - return np.array([]) - step = int(round(fs)) - step = max(step, 1) - return np.array([np.nanmean(env[i * step:(i + 1) * step]) - for i in range(nsec)]) - - a = _per_second(env_a, fs_a) - b = _per_second(env_b, fs_b) - n = min(len(a), len(b)) - if n < 3: - return dict(r=np.nan, n=n) - a = a[:n] - b = b[:n] - ok = np.isfinite(a) & np.isfinite(b) - if ok.sum() < 3 or a[ok].std() == 0 or b[ok].std() == 0: - return dict(r=np.nan, n=int(ok.sum())) - r = float(np.corrcoef(a[ok], b[ok])[0, 1]) - return dict(r=r, n=int(ok.sum())) +from micromotion.mocap import ( # noqa: F401 + compare_modality_envelopes, + dominant_frequency, + read_qtm_tsv, +) + +__all__ = [ + "compare_modality_envelopes", + "dominant_frequency", + "read_qtm_tsv", +] diff --git a/musicalgestures/_physio.py b/musicalgestures/_physio.py index 79115d4..c2eda23 100644 --- a/musicalgestures/_physio.py +++ b/musicalgestures/_physio.py @@ -1,141 +1,24 @@ -""" -Physiology signal features for standstill / micromotion studies. +"""Re-exported from the ``micromotion`` package. -Two pure numpy/scipy surfaces ported from the "still standing" study: +These functions used to live here. They were moved to ``micromotion`` on 2026-07-29 so that +one implementation of quantity of motion exists rather than two, and MGT now depends on that +package instead of carrying its own copy. Behaviour is unchanged: this module's tests pass +against ``micromotion`` unmodified. -* :func:`respiration_rate` -- windowed breathing rate (breaths per minute) - from a respiration waveform, via band-pass filtering and a Welch spectral - peak per window. -* :func:`spectral_band_fractions` -- the fraction of a signal's Welch power - falling in each of a set of caller-supplied named frequency bands. This is - the generic "cardiorespiratory QoM" spectral-composition diagnostic with - the heart-rate/respiration bands supplied by the caller, so the function - carries no dependency on any particular physiological sensor. +The dependency points this way round on purpose. ``micromotion`` needs only numpy, scipy and +pandas, so someone analysing accelerometer data does not have to install a computer-vision +stack; MGT already depends on ``ambiscape`` the same way, and neither of those packages +imports MGT. -Source: still standing study (Jensenius) -- Deichman / Equivital physiology -analyses. +Import from ``micromotion`` directly in new code. """ -import numpy as np - - -def respiration_rate(waveform, fs, *, band=(0.1, 0.6), window_s=30, - step_s=30): - """ - Windowed respiration rate (breaths per minute) from a breathing waveform. - - Each analysis window is band-pass filtered to the respiration band and - its dominant frequency is taken as the Welch spectral peak inside that - band; the rate is that frequency times 60. Windows advance by ``step_s`` - seconds. The default band ``(0.1, 0.6)`` Hz corresponds to about - 6-36 breaths/min. Each window must contain at least 15 seconds of valid - samples for spectral estimation. - - Source: still standing study (Jensenius), Deichman respiration analysis - (``compute_qom_resp``). - - Args: - waveform (np.ndarray): 1-D respiration/breathing waveform. - fs (float): Sampling rate in Hz. - band (tuple, optional): ``(low, high)`` respiration band in Hz. - Defaults to ``(0.1, 0.6)``. - window_s (float, optional): Window length in seconds. Defaults to 30. - step_s (float, optional): Hop between windows in seconds. Defaults to - 30. - - Returns: - dict: ``{"rate_bpm", "times_s", "median_bpm"}`` where ``rate_bpm`` is - the per-window rate (breaths/min, ``nan`` for windows without a - clear peak), ``times_s`` the window centre times in seconds, and - ``median_bpm`` the median across valid windows. - """ - from scipy.signal import butter, filtfilt, welch - - x = np.asarray(waveform, dtype=float) - x = x[np.isfinite(x)] - lo, hi = band - nyq = fs / 2.0 - if len(x) < int(fs * window_s): - # single short window: still attempt one estimate - window_s = max(1.0, len(x) / fs) - - b, a = butter(2, [lo / nyq, hi / nyq], btype="band") - xf = filtfilt(b, a, x - x.mean()) - - win = int(fs * window_s) - step = int(fs * step_s) - win = max(win, 1) - step = max(step, 1) - - rates = [] - times = [] - for start in range(0, max(len(xf) - win + 1, 1), step): - seg = xf[start:start + win] - if len(seg) < int(fs * 15): # need at least 15 s for spectral estimation - rates.append(np.nan) - times.append((start + win / 2) / fs) - continue - nperseg = min(len(seg), int(fs * window_s)) - f, P = welch(seg, fs, nperseg=nperseg) - mask = (f >= lo) & (f <= hi) - if mask.any() and P[mask].sum() > 0: - fpk = f[mask][np.argmax(P[mask])] - rates.append(float(fpk * 60.0)) - else: - rates.append(np.nan) - times.append((start + win / 2) / fs) - - rate_bpm = np.array(rates, dtype=float) - return dict(rate_bpm=rate_bpm, times_s=np.array(times, dtype=float), - median_bpm=float(np.nanmedian(rate_bpm)) - if np.isfinite(rate_bpm).any() else np.nan) - - -def spectral_band_fractions(signal, fs, bands, *, total_band=(0.1, 8.0), - nperseg_s=20): - """ - Fraction of a signal's power in each of a set of named frequency bands. - - Estimates the Welch power spectrum and, for each named band in ``bands``, - returns that band's summed power divided by the summed power in - ``total_band``. This is the generic spectral-composition diagnostic used - for the "cardiorespiratory QoM artifact" analysis (e.g. how much of a - chest-accelerometer QoM signal sits in a cardiac vs a respiration band), - with the bands supplied by the caller so there is no built-in dependence - on a heart-rate or respiration sensor. Power is bin-summed on the Welch - grid; the study source integrated with trapz, which yields nearly - identical results on the uniform frequency spacing of Welch. - - Source: still standing study (Jensenius), Deichman chest-QoM - cardiorespiratory spectral-composition analysis (``deichman_full``). - - Args: - signal (np.ndarray): 1-D input signal. - fs (float): Sampling rate in Hz. - bands (dict): Mapping of band name to ``(low, high)`` in Hz, e.g. - ``{"cardiac": (0.9, 1.3), "resp": (0.12, 0.5)}``. - total_band (tuple, optional): ``(low, high)`` reference band whose - power is the denominator. Defaults to ``(0.1, 8.0)``. - nperseg_s (float, optional): Welch segment length in seconds. - Defaults to 20. - - Returns: - dict: Mapping of each band name to its power fraction in ``[0, 1]`` - (``nan`` if the total band contains no power). - """ - from scipy.signal import welch +from micromotion.physio import ( # noqa: F401 + respiration_rate, + spectral_band_fractions, +) - x = np.asarray(signal, dtype=float) - x = x[np.isfinite(x)] - if len(x) < 8: - return {name: np.nan for name in bands} - nperseg = min(len(x), max(8, int(fs * nperseg_s))) - f, P = welch(x - x.mean(), fs, nperseg=nperseg) - tlo, thi = total_band - total = P[(f >= tlo) & (f < thi)].sum() - if total <= 0: - return {name: np.nan for name in bands} - out = {} - for name, (lo, hi) in bands.items(): - out[name] = float(P[(f >= lo) & (f < hi)].sum() / total) - return out +__all__ = [ + "respiration_rate", + "spectral_band_fractions", +] diff --git a/musicalgestures/_posture.py b/musicalgestures/_posture.py index 4751bc0..88ee353 100644 --- a/musicalgestures/_posture.py +++ b/musicalgestures/_posture.py @@ -1,633 +1,44 @@ -""" -Posturography and standstill-sway metrics for centre-of-pressure (CoP) and -head/marker position signals. - -This module ports the "still standing" study's posturography stack into -pure numpy/scipy surfaces that operate on plain arrays -- no study-specific -loaders, axis conventions, or marker loops. Three families of measures are -provided: +"""Re-exported from the ``micromotion`` package. -* **Sway amount / geometry** -- :func:`cop_sway_metrics`, - :func:`confidence_ellipse_area`, :func:`convex_hull_area`. -* **Control dynamics / complexity** -- :func:`stabilogram_diffusion` - (Collins-De Luca SDA), :func:`dfa` (detrended fluctuation analysis), - :func:`sample_entropy`, :func:`spectral_edges`, :func:`sway_texture`, - :func:`principal_axis_projection`. -* **Direction / extent** -- :func:`sway_orientation`, :func:`axial_rayleigh`, - :func:`spatial_extent`. +These functions used to live here. They were moved to ``micromotion`` on 2026-07-29 so that +one implementation of quantity of motion exists rather than two, and MGT now depends on that +package instead of carrying its own copy. Behaviour is unchanged: this module's tests pass +against ``micromotion`` unmodified. -The from-scratch SDA / DFA / sample-entropy implementations are validated in -the test-suite against known-answer synthetic signals (white noise -> -DFA alpha ~= 0.5 and SDA Hurst ~= 0.5; a sine -> low sample entropy relative -to its shuffle). +The dependency points this way round on purpose. ``micromotion`` needs only numpy, scipy and +pandas, so someone analysing accelerometer data does not have to install a computer-vision +stack; MGT already depends on ``ambiscape`` the same way, and neither of those packages +imports MGT. -Source: still standing study (Jensenius) -- posturography and micromotion -analyses of the international "standstill" championships and related datasets. +Import from ``micromotion`` directly in new code. """ -import numpy as np - - -# --------------------------------------------------------------------------- -# Sway amount / geometry -# --------------------------------------------------------------------------- -def confidence_ellipse_area(xy, conf=0.95): - """ - Area of the confidence ellipse of a 2-D point cloud (e.g. a - centre-of-pressure trace). - - The ellipse is the standard bivariate-Gaussian confidence region - ``area = pi * chi2_conf,2df * sqrt(det Cov)`` where ``Cov`` is the - 2x2 covariance of the (mean-removed) points. For a CoP sway path this - is the classic 95% "sway-ellipse area". - - Source: still standing study (Jensenius), HpSp balance analysis. - - Args: - xy (np.ndarray): Point cloud of shape ``(T, 2)``. - conf (float, optional): Confidence level in ``(0, 1)``. Defaults to - 0.95. - - Returns: - float: Ellipse area in squared position units (e.g. mm^2), or - ``nan`` if fewer than three finite points are available. - """ - from scipy.stats import chi2 - - xy = np.asarray(xy, dtype=float) - xy = xy[np.isfinite(xy).all(axis=1)] - if len(xy) < 3: - return np.nan - cov = np.cov(xy.T) - return float(np.pi * chi2.ppf(conf, 2) * np.sqrt(max(np.linalg.det(cov), 0.0))) - - -def convex_hull_area(xy): - """ - Area of the 2-D convex hull of a point cloud. - - A non-parametric alternative to :func:`confidence_ellipse_area` for the - region occupied by a sway path: it makes no Gaussian assumption and is - driven by the outermost excursions. - - Source: still standing study (Jensenius); complements the confidence - ellipse used in the balance reports. - - Args: - xy (np.ndarray): Point cloud of shape ``(T, 2)``. - - Returns: - float: Convex-hull area in squared position units, or ``nan`` if - fewer than three non-collinear finite points are available. - """ - from scipy.spatial import ConvexHull, QhullError - - xy = np.asarray(xy, dtype=float) - xy = xy[np.isfinite(xy).all(axis=1)] - if len(xy) < 3: - return np.nan - try: - return float(ConvexHull(xy).volume) # 2-D "volume" == area - except QhullError: - return np.nan - - -def cop_sway_metrics(xy, t=None, fs=None, *, freq_band=(0.1, 5.0), - resample_fs=50.0): - """ - Standard centre-of-pressure (CoP) sway metrics from a 2-D sway path. - - Computes the classic posturographic descriptors: CoP path length and - path rate, the 95% confidence-ellipse area, medio-lateral (ML) and - antero-posterior (AP) ranges and standard deviations, the AP/ML range - and SD ratios, and the mean sway frequency of each axis (the - power-weighted mean frequency of a Welch spectrum inside ``freq_band``, - computed on a uniform grid at ``resample_fs``). - - The first column of ``xy`` is treated as ML and the second as AP, - matching the study convention. Sampling time may be given either as an - explicit time vector ``t`` (seconds; may be irregular) or a constant - rate ``fs`` (Hz); if neither is supplied a rate of 1 Hz is assumed. - - Source: still standing study (Jensenius), HpSp balance analysis - (``analyze_balance``). - - Args: - xy (np.ndarray): CoP path of shape ``(T, 2)`` as ``[ML, AP]`` in - position units (e.g. mm). - t (np.ndarray, optional): Per-sample timestamps in seconds. May be - irregular. Defaults to None. - fs (float, optional): Constant sampling rate in Hz, used when ``t`` - is not given. Defaults to None (interpreted as 1 Hz). - freq_band (tuple, optional): ``(low, high)`` band in Hz for the mean - sway frequency. Defaults to ``(0.1, 5.0)``. - resample_fs (float, optional): Uniform rate in Hz onto which the - path is interpolated before the spectral estimate. Defaults to - 50.0. - - Returns: - dict: Metrics with keys ``n``, ``dur``, ``fs_mean``, ``path_len``, - ``path_rate``, ``area95``, ``ml_range``, ``ap_range``, - ``ml_sd``, ``ap_sd``, ``ap_ml_range_ratio``, - ``ap_ml_sd_ratio``, ``mf_ml``, ``mf_ap`` and ``mf_mean``. - """ - from scipy import signal as _sig - - xy = np.asarray(xy, dtype=float) - if xy.ndim != 2 or xy.shape[1] != 2: - raise ValueError("xy must have shape (T, 2)") - ml = xy[:, 0] - ap = xy[:, 1] - n = len(xy) - - if t is not None: - t = np.asarray(t, dtype=float) - t = t - t[0] - elif fs is not None: - t = np.arange(n) / float(fs) - else: - t = np.arange(n, dtype=float) - dur = float(t[-1]) if n > 1 else 0.0 - fs_mean = (n - 1) / dur if dur > 0 else np.nan - - ml_c = ml - np.nanmean(ml) - ap_c = ap - np.nanmean(ap) - - dpath = np.sqrt(np.diff(ml) ** 2 + np.diff(ap) ** 2) - path_len = float(np.nansum(dpath)) - path_rate = path_len / dur if dur > 0 else np.nan - - ml_range = float(np.nanmax(ml) - np.nanmin(ml)) - ap_range = float(np.nanmax(ap) - np.nanmin(ap)) - ml_sd = float(np.nanstd(ml, ddof=1)) - ap_sd = float(np.nanstd(ap, ddof=1)) - - area95 = confidence_ellipse_area(np.column_stack([ml_c, ap_c]), conf=0.95) - - lo, hi = freq_band - tu = np.arange(0, dur, 1.0 / resample_fs) if dur > 0 else np.array([]) - - def _mean_freq(x): - if len(tu) < 4: - return np.nan - xu = np.interp(tu, t, x) - xu = _sig.detrend(xu) - nperseg = min(len(xu), int(resample_fs * 20)) - if nperseg < 4: - return np.nan - f, P = _sig.welch(xu, fs=resample_fs, nperseg=nperseg) - band = (f >= lo) & (f <= hi) - if P[band].sum() <= 0: - return np.nan - return float(np.sum(f[band] * P[band]) / np.sum(P[band])) - - mf_ml = _mean_freq(ml_c) - mf_ap = _mean_freq(ap_c) - - return dict( - n=n, dur=dur, fs_mean=fs_mean, - path_len=path_len, path_rate=path_rate, area95=area95, - ml_range=ml_range, ap_range=ap_range, ml_sd=ml_sd, ap_sd=ap_sd, - ap_ml_range_ratio=ap_range / ml_range if ml_range else np.nan, - ap_ml_sd_ratio=ap_sd / ml_sd if ml_sd else np.nan, - mf_ml=mf_ml, mf_ap=mf_ap, mf_mean=float(np.nanmean([mf_ml, mf_ap])), - ) - - -# --------------------------------------------------------------------------- -# Control dynamics / complexity -# --------------------------------------------------------------------------- -def principal_axis_projection(xy): - """ - Project a 2-D (or N-D) point cloud onto its principal axis. - - Runs a PCA on the mean-removed points and returns the 1-D coordinate - along the direction of greatest variance -- the natural 1-D reduction of - a sway path used by the dynamics/complexity measures. The PCA eigenvector - sign is arbitrary; the projection may be globally flipped across calls or - datasets. - - Source: still standing study (Jensenius), sway-dynamics analysis. - - Args: - xy (np.ndarray): Point cloud of shape ``(T, D)`` (typically - ``D == 2``). - - Returns: - np.ndarray: 1-D projection of shape ``(T,)``. - """ - xy = np.asarray(xy, dtype=float) - xy = xy - xy.mean(axis=0) - cov = np.cov(xy.T) - w, v = np.linalg.eigh(cov) - return xy @ v[:, -1] # eigvecs ascending -> last is major axis - - -def stabilogram_diffusion(xy, fs, *, short_max_s=0.6, long_min_s=1.5, - n_lags=40): - """ - Collins-De Luca stabilogram-diffusion analysis (SDA) of a sway path. - - Fits the mean-square-displacement (MSD) curve - ``<[r(t+dt) - r(t)]^2>`` versus time-lag ``dt`` in log-log space and - reports a short-term and a long-term Hurst exponent (each ``slope / 2``) - plus the critical crossover time where the two regression lines - intersect. Persistent (open-loop) drift gives a short-term Hurst above - 0.5; anti-persistent (closed-loop correction) gives a long-term Hurst - below 0.5. - - Source: still standing study (Jensenius), sway-dynamics analysis; - method of Collins & De Luca (1993). - - Args: - xy (np.ndarray): Sway path of shape ``(T, D)`` (``D >= 1``). A 1-D - input of shape ``(T,)`` is accepted and treated as a single - axis. - fs (float): Sampling rate in Hz. - short_max_s (float, optional): Upper bound (s) of the short-term - fitting window. Defaults to 0.6. - long_min_s (float, optional): Lower bound (s) of the long-term - fitting window. Defaults to 1.5. - n_lags (int, optional): Number of log-spaced lags at which the MSD - is evaluated. Defaults to 40. - - Returns: - dict: ``{"H_short", "H_long", "crossover_s"}``. Entries are ``nan`` - when a window contains fewer than three usable lags. - """ - xy = np.asarray(xy, dtype=float) - if xy.ndim == 1: - xy = xy[:, None] - n = len(xy) - if n < 8: - return dict(H_short=np.nan, H_long=np.nan, crossover_s=np.nan) - - lags = np.unique( - np.round(np.logspace(0, np.log10(max(n // 4, 2)), n_lags)).astype(int)) - lags = lags[lags >= 1] - msd = np.array([np.mean(np.sum((xy[l:] - xy[:-l]) ** 2, axis=1)) - for l in lags]) - t = lags / fs - ok = msd > 0 - t, msd = t[ok], msd[ok] - short = t < short_max_s - long = t > long_min_s - - def _slope_intercept(mask): - if mask.sum() < 3: - return np.nan, np.nan - a, b = np.polyfit(np.log(t[mask]), np.log(msd[mask]), 1) - return a, b - - a_s, b_s = _slope_intercept(short) - a_l, b_l = _slope_intercept(long) - H_short = a_s / 2.0 if np.isfinite(a_s) else np.nan - H_long = a_l / 2.0 if np.isfinite(a_l) else np.nan - crossover = np.nan - if np.isfinite(a_s) and np.isfinite(a_l) and abs(a_s - a_l) > 1e-9: - crossover = float(np.exp((b_l - b_s) / (a_s - a_l))) - return dict(H_short=float(H_short), H_long=float(H_long), - crossover_s=crossover) - - -def dfa(x, *, n_scales=18, min_scale=10): - """ - Detrended fluctuation analysis (DFA) scaling exponent. - - Integrates the mean-removed signal, then measures the RMS of the - linearly-detrended integrated profile within non-overlapping windows of - increasing size; the slope of ``log F(n)`` versus ``log n`` is the DFA - exponent ``alpha``. White noise gives ``alpha ~= 0.5``; a random walk - (Brownian) gives ``alpha ~= 1.5``; ``alpha == 1`` is 1/f noise. - - Source: still standing study (Jensenius), sway-complexity analysis; - method of Peng et al. (1994). - - Args: - x (np.ndarray): 1-D input signal. - n_scales (int, optional): Number of log-spaced window sizes. - Defaults to 18. - min_scale (int, optional): Smallest window size in samples. Defaults - to 10. - - Returns: - float: The DFA exponent ``alpha`` (``nan`` if too short). - """ - x = np.asarray(x, dtype=float) - x = x[np.isfinite(x)] - n = len(x) - if n < 4 * min_scale: - return np.nan - y = np.cumsum(x - x.mean()) - scales = np.unique( - np.round(np.logspace(np.log10(min_scale), - np.log10(n // 4), n_scales)).astype(int)) - F = [] - used = [] - for m in scales: - k = n // m - if k < 1: - continue - seg = y[:k * m].reshape(k, m) - tt = np.arange(m) - rms = [np.sqrt(np.mean((s - np.polyval(np.polyfit(tt, s, 1), tt)) ** 2)) - for s in seg] - fm = np.mean(rms) - if fm > 0: - F.append(fm) - used.append(m) - if len(F) < 2: - return np.nan - return float(np.polyfit(np.log(used), np.log(F), 1)[0]) - - -def sample_entropy(x, m=2, r=0.2): - """ - Sample entropy (SampEn) of a 1-D signal. - - Measures regularity/predictability: the negative log conditional - probability that sequences close (within tolerance ``r``) for ``m`` - samples remain close for ``m + 1`` samples, self-matches excluded. - Lower values mean more repetitive/predictable signals. The signal is - z-scored internally so ``r`` is expressed as a fraction of its standard - deviation. Neighbour counts use a Chebyshev (max-norm) KD-tree. - - Source: still standing study (Jensenius), sway-complexity analysis; - method of Richman & Moorman (2000). - - Args: - x (np.ndarray): 1-D input signal. - m (int, optional): Embedding (template) length. Defaults to 2. - r (float, optional): Tolerance as a fraction of the signal SD. - Defaults to 0.2. - - Returns: - float: Sample entropy (``nan`` if undefined, e.g. no matches). - """ - from scipy.spatial import cKDTree - - x = np.asarray(x, dtype=float) - x = x[np.isfinite(x)] - N = len(x) - if N < m + 2: - return np.nan - x = (x - x.mean()) / (x.std() + 1e-12) - - def _count(mm): - emb = np.array([x[i:i + mm] for i in range(N - mm + 1)]) - tree = cKDTree(emb) - pairs = tree.query_ball_point(emb, r, p=np.inf) - return sum(len(p) - 1 for p in pairs) # exclude self-match - - B = _count(m) - A = _count(m + 1) - if B == 0 or A == 0: - return np.nan - return float(-np.log(A / B)) - - -def spectral_edges(x, fs, *, edges=(0.5, 0.95), nperseg=None): - """ - Spectral-edge frequencies of a signal. - - Returns the frequencies below which a given cumulative fraction of the - Welch power spectrum lies. With the default ``edges`` the first value is - the median frequency (50% edge) and the second the 95% spectral-edge - frequency, two standard descriptors of sway spectral shape. - - Source: still standing study (Jensenius), sway-dynamics analysis. - - Args: - x (np.ndarray): 1-D input signal. - fs (float): Sampling rate in Hz. - edges (tuple, optional): Cumulative-power fractions in ``(0, 1)``. - Defaults to ``(0.5, 0.95)``. - nperseg (int, optional): Welch segment length in samples. Defaults - to ``min(2048, len(x))``. - - Returns: - dict: Mapping ``"f"`` (e.g. ``"f50"``, ``"f95"``) to the edge - frequency in Hz. - """ - from scipy import signal as _sig - - x = np.asarray(x, dtype=float) - x = x[np.isfinite(x)] - if len(x) < 8: - return {f"f{int(round(e * 100))}": np.nan for e in edges} - if nperseg is None: - nperseg = min(2048, len(x)) - f, P = _sig.welch(x - x.mean(), fs, nperseg=nperseg) - cP = np.cumsum(P) - if cP[-1] <= 0: - return {f"f{int(round(e * 100))}": np.nan for e in edges} - cP = cP / cP[-1] - return {f"f{int(round(e * 100))}": float(np.interp(e, cP, f)) - for e in edges} - - -def sway_texture(speed, fs, *, frozen_threshold=2.0): - """ - Micro-texture of a sway speed signal: frozen fraction and burst rate. - - Distinguishes a smooth wander from intermittent ballistic corrections. - The frozen fraction is the share of time the speed is below - ``frozen_threshold``; the burst rate is the number of upward threshold - crossings (onset of a velocity burst) per minute. - - Source: still standing study (Jensenius), sway-texture analysis. - - Args: - speed (np.ndarray): 1-D speed signal (e.g. mm/s). - fs (float): Sampling rate in Hz. - frozen_threshold (float, optional): Speed below which the signal is - considered "frozen", in the units of ``speed``. Defaults to 2.0. - - Returns: - dict: ``{"frozen_fraction", "burst_rate"}`` where ``burst_rate`` is - in bursts per minute. - """ - speed = np.asarray(speed, dtype=float) - speed = speed[np.isfinite(speed)] - if len(speed) < 2: - return dict(frozen_fraction=np.nan, burst_rate=np.nan) - frozen = float((speed < frozen_threshold).mean()) - above = (speed >= frozen_threshold).astype(int) - bursts = int(np.sum(np.diff(above) == 1)) - minutes = len(speed) / fs / 60.0 - burst_rate = bursts / minutes if minutes > 0 else np.nan - return dict(frozen_fraction=frozen, burst_rate=burst_rate) - - -# --------------------------------------------------------------------------- -# Direction / extent -# --------------------------------------------------------------------------- -def sway_orientation(xy): - """ - Principal sway-axis orientation and anisotropy of a 2-D point cloud. - - A PCA of the mean-removed horizontal positions gives the orientation of - the major axis as an axial angle in ``[0, 180)`` degrees (undirected -- - a line, not an arrow) and the anisotropy ``sqrt(lambda_max / lambda_min)`` - of the sway ellipse. Anisotropy 1.0 is isotropic/circular; values above - ~1.3 indicate clearly directional sway. - - Source: still standing study (Jensenius), sway-direction analysis. - - Args: - xy (np.ndarray): Point cloud of shape ``(T, 2)``. - - Returns: - dict: ``{"angle_deg", "anisotropy"}``. Both are ``nan`` if the - covariance is degenerate or there are too few finite points. - """ - xy = np.asarray(xy, dtype=float) - xy = xy[np.isfinite(xy).all(axis=1)] - if len(xy) < 3: - return dict(angle_deg=np.nan, anisotropy=np.nan) - x = xy[:, 0] - xy[:, 0].mean() - y = xy[:, 1] - xy[:, 1].mean() - cov = np.cov(x, y) - w, v = np.linalg.eigh(cov) # ascending eigenvalues - if w[0] <= 0: - return dict(angle_deg=np.nan, anisotropy=np.nan) - angle = float(np.degrees(np.arctan2(v[1, 1], v[0, 1])) % 180.0) - anisotropy = float(np.sqrt(w[1] / w[0])) - return dict(angle_deg=angle, anisotropy=anisotropy) - - -def axial_rayleigh(angles_deg): - """ - Axial Rayleigh test for a preferred orientation among axial angles. - - Tests whether a sample of axial (undirected, ``[0, 180)`` deg) angles -- - e.g. per-session principal sway axes -- clusters around a common - orientation. Angles are doubled to map the axial circle onto the full - circle before computing the mean resultant length ``R`` and the Rayleigh - p-value (small ``p`` with large ``R`` means a shared preferred axis). - - Source: still standing study (Jensenius), sway-direction analysis. - - Args: - angles_deg (np.ndarray): Axial angles in degrees. - - Returns: - dict: ``{"R", "p", "mean_axis_deg", "n"}``. - """ - a = 2 * np.radians(np.asarray(angles_deg, dtype=float)) - a = a[np.isfinite(a)] - n = len(a) - if n < 2: - return dict(R=np.nan, p=np.nan, mean_axis_deg=np.nan, n=n) - C = np.mean(np.cos(a)) - S = np.mean(np.sin(a)) - R = float(np.hypot(C, S)) - Z = n * R * R - p = float(np.exp(-Z) * (1 + (2 * Z - Z * Z) / (4 * n))) - mean_axis = float(np.degrees(0.5 * np.arctan2(S, C)) % 180.0) - return dict(R=R, p=p, mean_axis_deg=mean_axis, n=n) - - -def spatial_extent(pos, fs, *, ellipse_conf=0.95, window_s=20.0, - vertical_axis=None): - """ - Spatial extent / occupied volume of a 3-D (or 2-D) position trace. - - Complements sway magnitude (QoM) by describing *how large a region* a - marker occupies. Reports the RMS dispersion radius about the session - centroid, the Gaussian confidence-ellipsoid volume and its cube-root - radius, the mean within-window dispersion (which removes slow drift), - and a drift ratio ``full_dispersion / within_window_dispersion`` (``> 1`` - when slow drift enlarges the occupied region over the session). When a - ``vertical_axis`` is given the drift is additionally split into - horizontal and vertical components. - - Source: still standing study (Jensenius), spatial-range analysis - (``session_metrics``). - - Args: - pos (np.ndarray): Position trace of shape ``(T, D)`` with ``D`` 2 or - 3, in position units (e.g. mm). - fs (float): Sampling rate in Hz. - ellipse_conf (float, optional): Confidence level for the ellipsoid - volume. Defaults to 0.95. - window_s (float, optional): Window length in seconds for the - within-window dispersion. Defaults to 20.0. - vertical_axis (int, optional): Index of the vertical axis (``0``, - ``1`` or ``2``); when given, drift is decomposed into horizontal - and vertical parts. Defaults to None. - - Returns: - dict: ``dispersion``, ``ellipsoid_volume``, ``ellipsoid_radius``, - ``within_window_dispersion``, ``drift_ratio`` and (when - ``vertical_axis`` is set) ``drift_horizontal`` and - ``drift_vertical``. Returns ``None`` if fewer than one window of - finite samples is available. - """ - from scipy.stats import chi2 - - pos = np.asarray(pos, dtype=float) - if pos.ndim != 2: - raise ValueError("pos must have shape (T, D)") - D = pos.shape[1] - mask = np.isfinite(pos).all(axis=1) - P = pos[mask] - w = int(fs * window_s) - if len(P) < max(w, D + 2): - return None - - centroid = P.mean(axis=0) - disp_vec = P - centroid - d = np.sqrt((disp_vec * disp_vec).sum(axis=1)) - dispersion = float(np.sqrt((d * d).mean())) - - cov = np.cov(P.T) - det = np.linalg.det(cov) - chi_d = chi2.ppf(ellipse_conf, D) - # volume of a D-dim Gaussian confidence ellipsoid - if D == 2: - vol = np.pi * chi_d * np.sqrt(max(det, 0.0)) - else: # D == 3 - vol = (4.0 / 3.0) * np.pi * (chi_d ** 1.5) * np.sqrt(max(det, 0.0)) - ellipsoid_volume = float(vol) - ellipsoid_radius = float(vol ** (1.0 / D)) - - def _within_window(Q): - cols = Q.shape[1] if Q.ndim == 2 else 1 - Q = Q.reshape(len(Q), cols) - wr = [] - for s in range(0, len(Q) - w + 1, w): - seg = Q[s:s + w] - cc = seg.mean(axis=0) - e = seg - cc - wr.append(np.sqrt((e * e).sum(axis=1).mean())) - return np.mean(wr) if wr else np.nan - - def _full_dispersion(Q): - cols = Q.shape[1] if Q.ndim == 2 else 1 - Q = Q.reshape(len(Q), cols) - cc = Q.mean(axis=0) - e = Q - cc - return np.sqrt((e * e).sum(axis=1).mean()) - - win_disp = float(_within_window(P)) - drift_ratio = (dispersion / win_disp - if win_disp and win_disp > 0 else np.nan) - - out = dict(dispersion=dispersion, ellipsoid_volume=ellipsoid_volume, - ellipsoid_radius=ellipsoid_radius, - within_window_dispersion=win_disp, drift_ratio=drift_ratio) - - if vertical_axis is not None: - hax = [k for k in range(D) if k != vertical_axis] - - def _comp_drift(cols): - Q = P[:, cols] - wm = _within_window(Q) - full = _full_dispersion(Q) - return full / wm if wm and wm > 0 else np.nan - - out["drift_horizontal"] = float(_comp_drift(hax)) - out["drift_vertical"] = float(_comp_drift([vertical_axis])) - - return out +from micromotion.balance import ( # noqa: F401 + axial_rayleigh, + confidence_ellipse_area, + convex_hull_area, + cop_sway_metrics, + dfa, + principal_axis_projection, + sample_entropy, + spatial_extent, + spectral_edges, + stabilogram_diffusion, + sway_orientation, + sway_texture, +) + +__all__ = [ + "axial_rayleigh", + "confidence_ellipse_area", + "convex_hull_area", + "cop_sway_metrics", + "dfa", + "principal_axis_projection", + "sample_entropy", + "spatial_extent", + "spectral_edges", + "stabilogram_diffusion", + "sway_orientation", + "sway_texture", +] diff --git a/musicalgestures/_qom.py b/musicalgestures/_qom.py index a2f29ba..3e7f2d1 100644 --- a/musicalgestures/_qom.py +++ b/musicalgestures/_qom.py @@ -1,384 +1,38 @@ -""" -Quantity-of-motion (QoM) signal cores for position, pose and accelerometer -data. +"""Re-exported from the ``micromotion`` package. -Band-limited QoM (with an automatic decimate+SOS regime for very low -frequency bands), accelerometer-to-speed integration, per-landmark-group -pose QoM, body-scale normalisation for framing-invariant comparisons, -spatial grid QoM, and small envelope/binning helpers. +These functions used to live here. They were moved to ``micromotion`` on 2026-07-29 so that +one implementation of quantity of motion exists rather than two, and MGT now depends on that +package instead of carrying its own copy. Behaviour is unchanged: this module's tests pass +against ``micromotion`` unmodified. -These functions are independent of the MgVideo/MgAudio classes and operate -on plain numpy arrays (marker/landmark trajectories, accelerometer data, -grayscale frame stacks, 1-D signals). +The dependency points this way round on purpose. ``micromotion`` needs only numpy, scipy and +pandas, so someone analysing accelerometer data does not have to install a computer-vision +stack; MGT already depends on ``ambiscape`` the same way, and neither of those packages +imports MGT. -Sources: stillstanding study and Westney-comparisons study (Jensenius). +Import from ``micromotion`` directly in new code. """ -import numpy as np - - -def _interp_nans(x): - """Linearly interpolate non-finite entries of a 1-D array (returns a copy).""" - x = np.asarray(x, float).copy() - m = np.isfinite(x) - if m.sum() > 2 and not m.all(): - x = np.interp(np.arange(len(x)), np.flatnonzero(m), x[m]) - return x - - -def envelope(x, fs, smooth=1.0, normalize=True): - """ - Smooth, optionally z-scored envelope of a signal: Savitzky-Golay - smoothing (order 2, window `smooth` seconds) followed by - standardisation. Used to compare motion/audio envelopes across sources - on a common, amplitude-free scale. - - Source: Westney-comparisons study (Jensenius). - - Args: - x (np.ndarray): Input 1-D signal. - fs (float): Sampling rate of the signal (Hz). - smooth (float, optional): Smoothing window in seconds. None or 0 disables - smoothing. Defaults to 1.0. - normalize (bool, optional): If True, z-score the result. Defaults to True. - - Returns: - np.ndarray: The smoothed (and optionally z-scored) envelope, same length - as the input. - """ - from scipy.signal import savgol_filter - x = np.asarray(x, float) - if smooth: - w = max(3, int(smooth * fs) | 1) - if len(x) > w: - x = savgol_filter(x, w, 2) - if normalize: - x = (x - x.mean()) / (x.std() + 1e-9) - return x - - -def bin_series(x, fs, bin_s=1.0): - """ - Mean of consecutive, non-overlapping bins of a signal (e.g. a per-second - quantity-of-motion envelope from a per-frame speed series). Trailing - samples that do not fill a whole bin are dropped. - - Source: stillstanding study (Jensenius); also used in the - Westney-comparisons study as a per-second envelope. - - Args: - x (np.ndarray): Input 1-D signal. - fs (float): Sampling rate of the signal (Hz). - bin_s (float, optional): Bin length in seconds. Defaults to 1.0. - - Returns: - np.ndarray: One mean value per bin (empty if the signal is shorter than - two bins). - """ - x = np.asarray(x, float) - step = max(1, int(round(fs * bin_s))) - n = len(x) // step - if n < 2: - return np.array([]) - return x[:n * step].reshape(n, step).mean(axis=1) - - -def band_limited_qom(pos, fs, lo=0.3, hi=15.0, order=4, auto_decimate=True): - """ - Band-limited quantity of motion from a position trajectory: the position - is band-pass filtered (zero phase) to `[lo, hi]` Hz and the QoM is the - per-frame speed, i.e. the Euclidean norm of the first difference times - the sampling rate (units of the input per second). - - For very low bands relative to the sampling rate (band edge below about - fs/40), a direct high-order band-pass is numerically fragile; in that - regime the trajectory is first decimated (zero phase) so the band sits - comfortably in the new Nyquist range, then filtered with a second-order - section (SOS) band-pass. This is the "slow sway" regime (e.g. 0.1-0.5 Hz - postural sway from 100 Hz mocap). Set `auto_decimate=False` to force the - direct filter. - - Source: stillstanding study and Westney-comparisons study (Jensenius) -- - this unifies the band-limited QoM cores used on mocap markers (mm), - MediaPipe landmarks (px) and slow postural sway across both studies. - - Args: - pos (np.ndarray): Position trajectory of shape (N,) or (N, D) (e.g. D=2 - image coordinates or D=3 mocap coordinates). Non-finite samples are - linearly interpolated per dimension. - fs (float): Sampling rate of the trajectory (Hz). - lo (float, optional): Lower band edge (Hz). Defaults to 0.3. - hi (float, optional): Upper band edge (Hz), clipped to 0.9 x Nyquist. - Defaults to 15.0. - order (int, optional): Butterworth order of the direct band-pass. - Defaults to 4. - auto_decimate (bool, optional): Enable the decimate+SOS low-band regime. - Defaults to True. - - Returns: - tuple: `(speed, fs_out)` where `speed` is the per-frame speed series - (length N-1, or shorter when decimated) and `fs_out` is its - sampling rate (equal to `fs` unless decimated). `speed` is empty - (and `fs_out` equals the input `fs`) when the input has fewer - than `int(fs) + 5` samples, or when it still contains non-finite - samples after per-dimension interpolation (i.e. a dimension had - fewer than 3 finite samples to interpolate from). In the - auto-decimate regime, `speed` is also empty (with `fs_out` the - decimated rate) when decimation leaves fewer than ~30 samples -- - too few for a stable SOS band-pass. - - Raises: - ValueError: If the band is invalid, i.e. does not satisfy - `0 < lo < hi <= 0.45*fs` (after `hi` is clipped to 0.9 x Nyquist). - """ - from scipy import signal - pos = np.asarray(pos, float) - if pos.ndim == 1: - pos = pos[:, None] - pos = np.column_stack([_interp_nans(pos[:, i]) for i in range(pos.shape[1])]) - - hi_eff = min(hi, 0.9 * fs / 2) - if not (0 < lo < hi_eff <= 0.45 * fs + 1e-9): - raise ValueError("band must satisfy 0 < lo < hi <= 0.45*fs") - - if len(pos) < int(fs) + 5 or not np.isfinite(pos).all(): - return np.array([]), fs - - if auto_decimate and fs / hi_eff >= 40: - q = int(min(13, fs // (20 * hi_eff))) - pos = signal.decimate(pos, q, axis=0, zero_phase=True) - fs_out = fs / q - if len(pos) < 30: - return np.array([]), fs_out - sos = signal.butter(2, [lo / (fs_out / 2), hi_eff / (fs_out / 2)], - btype="band", output="sos") - filtered = signal.sosfiltfilt(sos, pos, axis=0) - else: - fs_out = fs - b, a = signal.butter(order, [lo / (fs / 2), hi_eff / (fs / 2)], btype="band") - filtered = signal.filtfilt(b, a, pos, axis=0) - speed = np.linalg.norm(np.diff(filtered, axis=0), axis=1) * fs_out - return speed, fs_out - - -def accel_to_speed(acc, fs, highpass=0.3, order=2, normalize_gravity=False): - """ - Integrated speed from a 3-axis accelerometer: each axis is high-pass - filtered (removing gravity and DC), integrated to velocity, high-pass - filtered again (killing integration drift), and the speed is the - Euclidean norm of the velocity (m/s for input in m/s^2). - - Source: stillstanding study (Jensenius) -- the "corpus method" for - integrated quantity of motion from chest-worn accelerometers. - - Args: - acc (np.ndarray): Acceleration of shape (N, 3) in m/s^2 (or raw counts - with `normalize_gravity=True`). - fs (float): Sampling rate (Hz). - highpass (float, optional): High-pass cutoff (Hz) used both before and - after integration. Defaults to 0.3. - order (int, optional): Butterworth order of the high-pass filters. - Defaults to 2. - normalize_gravity (bool, optional): If True, rescale the raw input so - that the median vector magnitude equals 1 g (9.80665 m/s^2) before - filtering -- useful for uncalibrated sensors whose resting output - should be gravity. Defaults to False. - - Returns: - np.ndarray: Speed series of length N (m/s). - """ - from scipy.signal import butter, filtfilt - G = 9.80665 - acc = np.asarray(acc, float) - if normalize_gravity: - norm = np.linalg.norm(acc, axis=1) - acc = acc / np.median(norm) * G - b, a = butter(order, highpass / (fs / 2), btype="high") - acc = filtfilt(b, a, acc, axis=0) - vel = np.cumsum(acc, axis=0) / fs - vel = filtfilt(b, a, vel, axis=0) - return np.linalg.norm(vel, axis=1) - - -def group_qom(points, fs, lo=0.3, hi=15.0, **kwargs): - """ - Mean band-limited quantity of motion over a group of markers/landmarks, - plus the group's mean speed envelope: each trajectory is passed through - `band_limited_qom` and the per-trajectory speeds are averaged. - - Source: stillstanding study and Westney-comparisons study (Jensenius) -- - per-body-part QoM (head, shoulders, arms, wrists) from mocap markers and - pose landmarks. - - Args: - points (np.ndarray): Trajectories of shape (N, M, D): N frames, M - markers/landmarks, D spatial dimensions. - fs (float): Sampling rate (Hz). - lo (float, optional): Lower band edge (Hz). Defaults to 0.3. - hi (float, optional): Upper band edge (Hz). Defaults to 15.0. - **kwargs: Passed on to `band_limited_qom`. - - Returns: - tuple: `(qom, speed, fs_out)` where `qom` is the mean speed across - markers and time (NaN if no marker yields a valid series), `speed` - is the group's mean per-frame speed series, and `fs_out` its - sampling rate. - """ - points = np.asarray(points, float) - speeds, fs_out = [], fs - for m in range(points.shape[1]): - sp, fs_out = band_limited_qom(points[:, m, :], fs, lo=lo, hi=hi, **kwargs) - if len(sp): - speeds.append(sp) - if not speeds: - return np.nan, np.array([]), fs_out - L = min(len(s) for s in speeds) - mean_speed = np.mean([s[:L] for s in speeds], axis=0) - return float(np.mean([s.mean() for s in speeds])), mean_speed, fs_out - - -def pose_qom(landmarks, fs, lo=0.3, hi=5.0, **kwargs): - """ - Band-limited quantity of motion of 2-D pose landmarks (px/s): a thin - wrapper around `group_qom` with the band used for image-space pose - trajectories (0.3-5 Hz), where higher bands are dominated by landmark - jitter rather than motion. - - Source: Westney-comparisons study (Jensenius). - - Args: - landmarks (np.ndarray): Landmark trajectories of shape (N, L, 2) in - pixels (a single landmark of shape (N, 2) is also accepted). - fs (float): Sampling rate (Hz, e.g. video frame rate). - lo (float, optional): Lower band edge (Hz). Defaults to 0.3. - hi (float, optional): Upper band edge (Hz). Defaults to 5.0. - **kwargs: Passed on to `band_limited_qom`. - - Returns: - tuple: `(qom, speed, fs_out)` as in `group_qom`. - """ - landmarks = np.asarray(landmarks, float) - if landmarks.ndim == 2: - landmarks = landmarks[:, None, :] - return group_qom(landmarks, fs, lo=lo, hi=hi, **kwargs) - - -def body_scale(landmarks, upper=(11, 12), lower=(23, 24)): - """ - Body-size scale (in the landmarks' own units, e.g. pixels) as the median - torso length: the distance from the midpoint of the `upper` landmarks - (shoulders) to the midpoint of the `lower` landmarks (hips). The torso - length is preferred over shoulder width because it stays robust in a - profile view, where the shoulder width collapses. - - The default indices are MediaPipe Pose landmarks (11/12 shoulders, - 23/24 hips). - - Source: Westney-comparisons study (Jensenius). - - Args: - landmarks (np.ndarray): Landmark trajectories of shape (N, L, C) with - C >= 2; only the first two coordinates are used. - upper (tuple, optional): Indices of the two shoulder landmarks. - Defaults to (11, 12). - lower (tuple, optional): Indices of the two hip landmarks. - Defaults to (23, 24). - - Returns: - float: Median torso length (NaN if no finite frames). - """ - landmarks = np.asarray(landmarks, float) - um = (landmarks[:, upper[0], :2] + landmarks[:, upper[1], :2]) / 2 - lm = (landmarks[:, lower[0], :2] + landmarks[:, lower[1], :2]) / 2 - d = np.linalg.norm(um - lm, axis=1) - d = d[np.isfinite(d)] - return float(np.median(d)) if len(d) else np.nan - - -def normalized_qom(landmarks, fs, scale=None, lo=0.3, hi=5.0, - upper=(11, 12), lower=(23, 24), **kwargs): - """ - Body-scale-normalised quantity of motion (body-lengths per second): - the pose QoM divided by the performer's own body scale (median torso - length, see `body_scale`). Being dimensionless, this is invariant to - camera framing/zoom and comparable across recordings. - - Source: Westney-comparisons study (Jensenius) -- framing-invariant - with/without-audience comparison of a pianist's motion. - - Args: - landmarks (np.ndarray): Landmark trajectories of shape (N, L, 2). - fs (float): Sampling rate (Hz). - scale (float, optional): Precomputed body scale. Defaults to None (which - computes `body_scale(landmarks, upper, lower)`). - lo (float, optional): Lower band edge (Hz). Defaults to 0.3. - hi (float, optional): Upper band edge (Hz). Defaults to 5.0. - upper (tuple, optional): Shoulder landmark indices for `body_scale`. Defaults to (11, 12). - lower (tuple, optional): Hip landmark indices for `body_scale`. Defaults to (23, 24). - **kwargs: Passed on to `band_limited_qom`. - - Returns: - tuple: `(qom, speed, fs_out)` as in `group_qom`, with both `qom` and - `speed` divided by the body scale. When `scale` is non-finite - (e.g. `body_scale` found no finite torso-length sample) or not - strictly positive (degenerate, coincident upper/lower landmarks), - division would otherwise silently propagate NaN/inf through - `qom` and `speed`; instead both are explicitly returned as NaN - (`qom` as a NaN scalar, `speed` as an all-NaN array of the same - shape) so the invalid-scale case is unambiguous rather than - merely inferred from the arithmetic. - """ - landmarks = np.asarray(landmarks, float) - if scale is None: - scale = body_scale(landmarks, upper=upper, lower=lower) - qom, speed, fs_out = pose_qom(landmarks, fs, lo=lo, hi=hi, **kwargs) - if not np.isfinite(scale) or scale <= 0: - return float("nan"), np.full_like(speed, np.nan), fs_out - return qom / scale, speed / scale, fs_out - - -def grid_qom(frames, grid=(6, 4), region=(0.0, 1.0, 0.0, 1.0), threshold=8.0): - """ - Spatial grid quantity of motion from a stack of grayscale frames: the - absolute inter-frame difference is thresholded (small differences set to - zero to suppress sensor noise) and averaged within each cell of a - `grid[0]` x `grid[1]` grid laid over `region`, yielding one motion time - series per cell plus a per-cell mean-motion heatmap. - - Source: Westney-comparisons study (Jensenius) -- audience-region motion - mapping in a concert hall. - - Args: - frames (np.ndarray): Grayscale frames of shape (T, H, W). - grid (tuple, optional): Grid size (columns, rows). Defaults to (6, 4). - region (tuple, optional): Region of interest as fractions - (x0, x1, y0, y1) of the frame. Defaults to the full frame. - threshold (float, optional): Absolute-difference threshold below which - pixel changes are zeroed (0-255 scale). Defaults to 8.0. - - Returns: - tuple: `(series, heat)` where `series` has shape (T-1, rows*cols) - (cells in row-major order) and `heat` has shape (rows, cols) with - each cell's time-mean motion. - - Raises: - ValueError: If `frames` is not 3-D (T, H, W), as in - `_motionanalysis.motiongram_data`. - """ - frames = np.asarray(frames, dtype=np.float32) - if frames.ndim != 3: - raise ValueError("grid_qom expects frames of shape (T, H, W)") - T, H, W = frames.shape - gx, gy = grid - x0, x1, y0, y1 = region - xs = np.linspace(int(x0 * W), int(x1 * W), gx + 1).astype(int) - ys = np.linspace(int(y0 * H), int(y1 * H), gy + 1).astype(int) - d = np.abs(np.diff(frames, axis=0)) - d[d < threshold] = 0.0 - series = np.empty((T - 1, gy * gx), dtype=np.float32) - for r in range(gy): - for c in range(gx): - cell = d[:, ys[r]:ys[r + 1], xs[c]:xs[c + 1]] - series[:, r * gx + c] = cell.mean(axis=(1, 2)) - heat = series.mean(axis=0).reshape(gy, gx) - return series, heat +from micromotion.qom import ( # noqa: F401 + accel_to_speed, + band_limited_qom, + bin_series, + body_scale, + envelope, + grid_qom, + group_qom, + normalized_qom, + pose_qom, +) + +__all__ = [ + "accel_to_speed", + "band_limited_qom", + "bin_series", + "body_scale", + "envelope", + "grid_qom", + "group_qom", + "normalized_qom", + "pose_qom", +] diff --git a/pyproject.toml b/pyproject.toml index e70abbb..f9761f9 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -24,6 +24,7 @@ classifiers = [ "Programming Language :: Python :: 3.12", ] dependencies = [ + "micromotion>=0.3", "numpy", "pandas", "matplotlib",