feat: implement compute_path_sinuosity core logic (draft)

movement/kinematics/__init__.py

@@ -18,6 +18,7 @@
 from movement.kinematics.path import (
     compute_directional_change,
     compute_path_length,
+    compute_path_sinuosity,
     compute_path_straightness,
     compute_turning_angle,
 )
@@ -31,6 +32,7 @@
     "compute_speed",
     "compute_path_length",
     "compute_path_straightness",
+    "compute_path_sinuosity",
     "compute_directional_change",
     "compute_time_derivative",
     "compute_pairwise_distances",

movement/kinematics/path.py

@@ -12,7 +12,10 @@
 import numpy as np
 import xarray as xr
 
-from movement.kinematics.kinematics import compute_backward_displacement
+from movement.kinematics.kinematics import (
+    compute_backward_displacement,
+    compute_forward_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
@@ -393,9 +396,173 @@ def compute_directional_change(
     return dc
 
 
+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 (Benhamou 2004).
+
+    Sinuosity quantifies the tortuosity of a random search 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 (Benhamou 2004, Eq. 8) 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 forward displacement vectors
+    via :func:`~movement.utils.vector.compute_norm` and
+    :func:`~movement.kinematics.compute_forward_displacement`.
+    Turning angles are computed via :func:`compute_turning_angle`,
+    which internally composes
+    :func:`~movement.kinematics.compute_backward_displacement` and
+    :func:`~movement.utils.vector.compute_signed_angle_2d`.
+
+    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``
+      via :func:`xarray.where` before the division is evaluated.
+    - :math:`\bar{p} = 0` (entirely stationary track): S is NaN,
+      propagated from the safe division guard on :math:`\bar{p}`.
+    - 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 = compute_norm(compute_forward_displacement(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 = np.cos(theta).mean(dim="time", skipna=True)  # type: ignore
+
+    # 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"
+    result.attrs["units"] = "1/sqrt(length)"
+
+    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.
 
@@ -405,6 +572,8 @@ def _validate_time_points(
         Position data with ``time`` and ``space`` dimensions.
     metric_name : str
         Used in the error message when there are fewer than 2 time points.
+    min_points : int, optional
+        The minimum number of time points required. Defaults to 2.
 
     Returns
     -------
@@ -414,7 +583,7 @@ 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 "