feat: implement compute_path_sinuosity core logic (draft)

movement/kinematics/__init__.py

@@ -19,6 +19,7 @@
     compute_directional_change,
     compute_path_deviation,
     compute_path_length,
+    compute_path_sinuosity,
     compute_path_straightness,
     compute_turning_angle,
 )
@@ -32,6 +33,7 @@
     "compute_speed",
     "compute_path_length",
     "compute_path_straightness",
+    "compute_path_sinuosity",
     "compute_directional_change",
     "compute_path_deviation",
     "compute_time_derivative",

movement/kinematics/path.py

@@ -12,7 +12,9 @@
 import numpy as np
 import xarray as xr
 
-from movement.kinematics.kinematics import compute_backward_displacement
+from movement.kinematics.kinematics import (
+    compute_backward_displacement,
+)
 from movement.utils.logging import logger
 from movement.utils.reports import report_nan_values
 from movement.utils.vector import compute_norm, compute_signed_angle_2d
@@ -510,18 +512,181 @@ def compute_path_deviation(
     return deviation
 
 
+def compute_path_sinuosity(
+    data: xr.DataArray,
+    nan_warn_threshold: float = 0.2,
+    min_step_length: float = 0.0,
+) -> xr.DataArray:
+    r"""Compute the sinuosity of a path.
+
+    Sinuosity (S) quantifies the tortuosity of a path by combining
+    turning angle statistics with step-length variability. Higher
+    values indicate more tortuous movement. A perfectly straight
+    path has S = 0.
+    ``
+
+    The corrected sinuosity index (Eq. 8 in [1]_) is defined as:
+
+    .. math::
+
+        S = 2\left[\bar{p}\left(
+            \frac{1+\bar{c}}{1-\bar{c}} + b^{2}
+        \right)\right]^{-1/2}
+
+    where :math:`\bar{p}` is the mean step length,
+    :math:`\bar{c} = \tfrac{1}{n}\sum_{i=1}^{n}\cos(\phi_i)` is the mean
+    cosine of turning angles, and
+    :math:`b = \mathrm{SD}(p_i)\,/\,\bar{p}` is the coefficient of
+    variation of step length.
+
+    Parameters
+    ----------
+    data : xarray.DataArray
+        The input data containing position information, with ``time``
+        and ``space`` (in Cartesian coordinates) as required dimensions.
+        To compute sinuosity over a specific time window, pre-slice with
+        ``data.sel(time=slice(start, stop))`` before passing in.
+    nan_warn_threshold : float, optional
+        If any point track in the data has at least (:math:`\ge`)
+        this proportion of values missing, a warning will be emitted.
+        Defaults to ``0.2`` (20%).
+    min_step_length : float, optional
+        Minimum step length threshold. Steps shorter than or equal to
+        this value are treated as stationary and excluded from all
+        summary statistics (:math:`\bar{p}`, :math:`\bar{c}`, :math:`b`).
+        Applied symmetrically here and passed to
+        :func:`compute_turning_angle`. Defaults to ``0.0``.
+
+    Returns
+    -------
+    xarray.DataArray
+        An xarray DataArray containing the computed sinuosity,
+        with dimensions matching those of the input data,
+        except ``time`` and ``space`` are removed.
+
+    See Also
+    --------
+    compute_path_length : Total path length between two time points.
+    compute_path_straightness : Net displacement divided by path length.
+    compute_turning_angle : Step-wise turning angle along a path.
+
+    Notes
+    -----
+    Step lengths are computed as the norm of backward displacement vectors
+    via :func:`~movement.utils.vector.compute_norm` and
+    :func:`~movement.kinematics.compute_backward_displacement`.
+    Turning angles are computed via :func:`compute_turning_angle`.
+
+    Steps shorter than or equal to ``min_step_length`` are masked to NaN
+    before computing any statistics. This mirrors the masking applied
+    inside :func:`compute_turning_angle`, ensuring noise steps are excluded
+    symmetrically from :math:`\bar{p}`, :math:`\bar{c}`, and :math:`b`.
+
+    NaN positions propagate to NaN step lengths and turning angles;
+    the statistics are then computed over the remaining valid samples
+    via ``skipna=True``.
+
+    Sinuosity has units of :math:`1/\sqrt{\text{length}}`, so its
+    numerical value depends on the position units of the input data.
+    Values are not directly comparable across datasets recorded in
+    different spatial units.
+
+    **Edge cases**:
+
+    - :math:`\bar{c} = 1` (perfectly straight path): the denominator
+      :math:`(1-\bar{c})` is zero; S is explicitly set to ``0.0``
+      before the division is evaluated.
+    - :math:`\bar{p} = 0` (entirely stationary track): S is NaN.
+    - All steps NaN: returns NaN.
+
+    References
+    ----------
+    .. [1] Benhamou, S. (2004). How to reliably estimate the tortuosity
+       of an animal's path: straightness, sinuosity, or fractal dimension?
+       *Journal of Theoretical Biology*, 229(2), 209–220.
+       https://doi.org/10.1016/j.jtbi.2004.03.016
+
+    Examples
+    --------
+    Compute sinuosity for all individuals and keypoints:
+
+    >>> sinuosity = compute_path_sinuosity(ds.position)
+
+    Constrain to a specific time window using xarray's ``.sel()``:
+
+    >>> sinuosity = compute_path_sinuosity(
+    ...     ds.position.sel(time=slice(10.0, 60.0))
+    ... )
+
+    Filter out sub-pixel noise steps:
+
+    >>> sinuosity = compute_path_sinuosity(ds.position, min_step_length=0.5)
+
+    """
+    data = _validate_time_points(
+        data, metric_name="path sinuosity", min_points=3
+    )
+
+    _warn_about_nan_proportion(data, nan_warn_threshold)
+
+    # --- Step lengths -------------------------------------------------------
+    step_lengths = _segment_lengths(data)
+
+    # Mask steps <= min_step_length symmetrically — matches the masking
+    # applied inside compute_turning_angle so noise steps are excluded
+    # from p_bar, c_bar, and b consistently.
+    step_lengths = step_lengths.where(step_lengths > min_step_length)
+
+    # --- Turning angles -----------------------------------------------------
+    theta = compute_turning_angle(data, min_step_length=min_step_length)
+
+    # --- Summary statistics (NaN-aware) -------------------------------------
+    p_bar = step_lengths.mean(dim="time", skipna=True)
+    c_bar = xr.apply_ufunc(np.cos, theta).mean(dim="time", skipna=True)
+
+    # Guard p_bar == 0 (entirely stationary track) before dividing to prevent
+    # RuntimeWarning. xarray-safe comparison preserves dimension labels.
+    is_stationary = abs(p_bar) < 1e-8
+    p_bar_safe = xr.where(is_stationary, np.nan, p_bar)
+
+    b = step_lengths.std(dim="time", skipna=True) / p_bar_safe
+
+    # --- Benhamou 2004 Eq. 8 ------------------------------------------------
+    # Guard c_bar == 1 before (1 - c_bar) is evaluated to avoid zero-division.
+    # xarray-safe: abs() keeps DataArray labels; np.isclose would drop them.
+    is_straight = abs(c_bar - 1.0) < 1e-8
+
+    angle_term = xr.where(
+        is_straight,
+        0.0,
+        (1.0 + c_bar) / (1.0 - c_bar),
+    )
+    result = 2.0 * (p_bar_safe * (angle_term + b**2)) ** -0.5
+
+    result = xr.where(is_straight, 0.0, result)
+
+    result.name = "sinuosity"
+    result.attrs["long_name"] = "Path Sinuosity"
+
+    return result
+
+
 def _validate_time_points(
     data: xr.DataArray,
     metric_name: str,
+    min_points: int = 2,
 ) -> xr.DataArray:
-    """Validate dims/coords and require at least 2 time points.
+    """Validate dims/coords and require at least ``min_points`` time points.
 
     Parameters
     ----------
     data : xarray.DataArray
         Position data with ``time`` and ``space`` dimensions.
     metric_name : str
-        Used in the error message when there are fewer than 2 time points.
+        Used in the error message when there are fewer than ``min_points``
+        time points.
+    min_points : int, optional
+        The minimum number of time points required. Defaults to 2.
 
     Returns
     -------
@@ -531,10 +696,10 @@ def _validate_time_points(
     """
     validate_dims_coords(data, {"time": [], "space": []})
     n_time = data.sizes["time"]
-    if n_time < 2:
+    if n_time < min_points:
         raise logger.error(
             ValueError(
-                "At least 2 time points are required to compute "
+                f"At least {min_points} time points are required to compute "
                 f"{metric_name}, but {n_time} were found."
             )
         )

tests/test_unit/test_kinematics/test_path.py

@@ -9,6 +9,7 @@
     compute_directional_change,
     compute_path_deviation,
     compute_path_length,
+    compute_path_sinuosity,
     compute_path_straightness,
     compute_turning_angle,
 )
@@ -801,3 +802,112 @@ def test_path_deviation_partially_degenerate_warns(straight_paths):
         assert np.isnan(result.sel(individual="id_0").values).all()
         id_1 = result.sel(individual="id_1")
         xr.testing.assert_allclose(id_1, xr.zeros_like(id_1))
+
+
+# ─────────────────────────────────────────────
+# Path sinuosity tests
+# ─────────────────────────────────────────────
+
+
+def test_path_sinuosity_too_few_timepoints(valid_poses_dataset):
+    """Test that sinuosity correctly enforces the minimum length requirement."""
+    # Slice the data to exactly 2 time points (which is only 1 segment)
+    position = valid_poses_dataset.position.isel(time=slice(0, 2))
+
+    with pytest.raises(ValueError):
+        compute_path_sinuosity(position)
+
+
+@pytest.mark.parametrize(
+    "fixture_name, expected_value",
+    [
+        pytest.param("straight_paths", 0.0, id="straight-line"),
+        pytest.param("stationary_paths", np.nan, id="stationary"),
+    ],
+)
+def test_path_sinuosity_known_values(request, fixture_name, expected_value):
+    """Test that sinuosity matches expected values for standard geometries."""
+    position = request.getfixturevalue(fixture_name)
+    result = compute_path_sinuosity(position)
+
+    if np.isnan(expected_value):
+        assert result.isnull().all()
+    else:
+        xr.testing.assert_allclose(
+            result,
+            xr.full_like(result, expected_value),
+            atol=1e-7,
+        )
+
+
+def test_path_sinuosity_perfect_reversals():
+    """Biological edge case: Trapped animal pacing.
+
+    A perfectly oscillating path (A -> B -> A -> B) with a constant stride.
+    Mean turning angle cosine (c̄) is -1, and step length variation (b) is 0.
+    Because tortuosity requires variance, sinuosity is mathematically undefined
+    here and should safely return NaN rather than evaluating to infinity.
+    """
+    xy = np.zeros((20, 2))
+    xy[1::2, 0] = 1.0  # Alternates x position exactly between 0 and 1
+    data = xr.DataArray(
+        xy,
+        dims=["time", "space"],
+        coords={"time": np.arange(20), "space": ["x", "y"]},
+    )
+    result = compute_path_sinuosity(data)
+
+    # Asserting it should safely return NaN without raising a division error
+    assert result.isnull().all()
+
+
+def test_path_sinuosity_variable_reversals():
+    """Biological edge case: Natural tortuous pacing.
+
+    An animal repeatedly reversing direction but with natural, variable
+    step lengths. This is genuinely tortuous movement and should compute
+    successfully to a finite, positive value.
+    """
+    rng = np.random.default_rng(11)
+    n = 30
+    x = np.zeros(n)
+    for i in range(1, n):
+        # Reverse direction each step, but multiply by a random stride length
+        x[i] = x[i - 1] + rng.choice([-1, 1]) * rng.uniform(0.5, 2.0)
+    xy = np.column_stack([x, np.zeros(n)])
+
+    data = xr.DataArray(
+        xy,
+        dims=["time", "space"],
+        coords={"time": np.arange(n), "space": ["x", "y"]},
+    )
+    result = compute_path_sinuosity(data)
+
+    assert not result.isnull().any()
+    assert float(result.item()) > 0
+    assert np.isfinite(result.item())
+
+
+def test_path_sinuosity_outlier_step():
+    """Artefact edge case: Tracking teleportation.
+
+    Simulates a tracker losing ID and assigning a point across the arena.
+    The massive jump creates an extreme standard deviation in step length (b).
+    The metric should not crash, but rather return a valid, finite float.
+    """
+    rng = np.random.default_rng(20)
+    n = 50
+    # Normal erratic movement
+    xy = np.column_stack([np.arange(n, dtype=float), rng.normal(0, 0.3, n)])
+    # Single massive tracking failure
+    xy[25] = [1000.0, 1000.0]
+
+    data = xr.DataArray(
+        xy,
+        dims=["time", "space"],
+        coords={"time": np.arange(n), "space": ["x", "y"]},
+    )
+    result = compute_path_sinuosity(data)
+
+    assert not result.isnull().any()
+    assert np.isfinite(result.item())