Extend smoothness functionality as raised in #66

Author: sanketsabharwal

pyeyesweb/low_level/smoothness.py

@@ -1,11 +1,11 @@
 """Movement smoothness analysis module.
 
-This module provides tools for quantifying the smoothness of movement velocity signals
+This module provides tools for quantifying the smoothness of movement signals
 using multiple metrics including SPARC (Spectral Arc Length) and Jerk RMS.
 Designed for real-time analysis of motion capture or sensor data.
 
-IMPORTANT: This module expects velocity signals as input, not position signals.
-Velocity should be computed from position data before analysis (for example velocity = sqrt(dx^2 + dy^2) / dt).
+The module can accept either position or velocity signals as input. When position
+is provided, velocity is automatically computed for SPARC analysis.
 
 Smoothness metrics are important indicators of movement quality in:
 1. Motor control assessment
@@ -18,17 +18,19 @@
 
 from pyeyesweb.data_models.sliding_window import SlidingWindow
 from pyeyesweb.utils.signal_processing import apply_savgol_filter
-from pyeyesweb.utils.math_utils import compute_sparc, compute_jerk_rms, normalize_signal
+from pyeyesweb.utils.math_utils import (
+    compute_sparc, compute_jerk_rms, normalize_signal, extract_velocity_from_position
+)
 from pyeyesweb.utils.validators import validate_numeric, validate_boolean
 
 
 class Smoothness:
-    """Compute movement smoothness metrics from velocity signal data.
+    """Compute movement smoothness metrics from position or velocity signal data.
 
     This class analyzes movement smoothness using SPARC (Spectral Arc Length)
-    and Jerk RMS metrics. This class expects velocity signals as input,
-    not position signals. It can optionally apply Savitzky-Golay filtering
-    to reduce noise before analysis.
+    and Jerk RMS metrics. It can accept either position or velocity signals as input.
+    When position is provided, velocity is automatically computed. It can optionally
+    apply Savitzky-Golay filtering to reduce noise before analysis.
 
     SPARC implementation is based on Balasubramanian et al. (2015) "On the analysis
     of movement smoothness" from Journal of NeuroEngineering and Rehabilitation.
@@ -44,26 +46,31 @@ class Smoothness:
         Sampling rate of the signal in Hz (default: 50.0).
     use_filter : bool, optional
         Whether to apply Savitzky-Golay filtering before analysis (default: True).
+    signal_type : {'velocity', 'position'}, optional
+        Type of input signal (default: 'velocity').
+        - 'velocity': Direct velocity input
+        - 'position': Position input (x, y coordinates or 1D position)
 
     Attributes
     ----------
     rate_hz : float
         Signal sampling rate.
     use_filter : bool
         Filter application flag.
+    signal_type : str
+        Type of input signal ('velocity' or 'position').
 
     Examples
     --------
     >>> from pyeyesweb.low_level.smoothness import Smoothness
     >>> from pyeyesweb.data_models.sliding_window import SlidingWindow
     >>> import numpy as np
     >>>
-    >>> # Generate sample velocity data
-    >>> # For example, from tracking hand movement: velocity = sqrt(dx^2 + dy^2) / dt
+    >>> # Example 1: Using velocity data (traditional approach)
     >>> t = np.linspace(0, 2, 200)
     >>> velocity_data = np.sin(2 * np.pi * t) + 0.1 * np.random.randn(200)
     >>>
-    >>> smooth = Smoothness(rate_hz=100.0, use_filter=True)
+    >>> smooth = Smoothness(rate_hz=100.0, use_filter=True, signal_type='velocity')
     >>> window = SlidingWindow(max_length=200, n_columns=1)
     >>>
     >>> # Add velocity data to the window
@@ -72,6 +79,22 @@ class Smoothness:
     >>>
     >>> result = smooth(window)
     >>> print(f"SPARC: {result['sparc']:.3f}, Jerk RMS: {result['jerk_rms']:.3f}")
+    >>>
+    >>> # Example 2: Using position data
+    >>> # 2D position data (x, y coordinates)
+    >>> t = np.linspace(0, 2, 200)
+    >>> x_positions = 10 * np.sin(2 * np.pi * t)
+    >>> y_positions = 5 * np.cos(2 * np.pi * t)
+    >>>
+    >>> smooth_pos = Smoothness(rate_hz=100.0, use_filter=True, signal_type='position')
+    >>> window_pos = SlidingWindow(max_length=200, n_columns=2)
+    >>>
+    >>> # Add position data to the window
+    >>> for x, y in zip(x_positions, y_positions):
+    ...     window_pos.append([x, y])
+    >>>
+    >>> result = smooth_pos(window_pos)
+    >>> print(f"SPARC: {result['sparc']:.3f}, Jerk RMS: {result['jerk_rms']:.3f}")
 
     Notes
     -----
@@ -86,10 +109,15 @@ class Smoothness:
     12(1), 1-11.
     """
 
-    def __init__(self, rate_hz=50.0, use_filter=True):
+    def __init__(self, rate_hz=50.0, use_filter=True, signal_type='velocity'):
         self.rate_hz = validate_numeric(rate_hz, 'rate_hz', min_val=0.01, max_val=100000)
         self.use_filter = validate_boolean(use_filter, 'use_filter')
 
+        # Validate signal_type
+        if signal_type not in ['velocity', 'position']:
+            raise ValueError(f"signal_type must be 'velocity' or 'position', got '{signal_type}'")
+        self.signal_type = signal_type
+
     def _filter_signal(self, signal):
         """Apply Savitzky-Golay filter if enabled and enough data.
 
@@ -108,37 +136,52 @@ def _filter_signal(self, signal):
         return apply_savgol_filter(signal, self.rate_hz)
 
     def __call__(self, sliding_window: SlidingWindow):
-        """Compute smoothness metrics from windowed velocity signal data.
+        """Compute smoothness metrics from windowed signal data.
 
         Parameters
         ----------
         sliding_window : SlidingWindow
-            Buffer containing velocity signal data to analyze.
-            Input should be velocity values, not position values.
+            Buffer containing signal data to analyze.
+            - If signal_type='velocity': expects velocity values
+            - If signal_type='position': expects position values (1D or 2D)
 
         Returns
         -------
         dict
             Dictionary containing:
             - 'sparc': Spectral Arc Length (closer to 0 = smoother).
                       Returns NaN if insufficient data.
-            - 'jerk_rms': RMS of jerk (computed from velocity input).
+            - 'jerk_rms': RMS of jerk.
                          Returns NaN if insufficient data.
         """
         if len(sliding_window) < 5:
             return {"sparc": np.nan, "jerk_rms": np.nan}
 
         signal, _ = sliding_window.to_array()
 
-        # If multi-channel, compute for first channel only
-        if signal.ndim > 1 and signal.shape[1] > 1:
-            signal = signal[:, 0]
-
-        filtered = self._filter_signal(signal.squeeze())
-        normalized = normalize_signal(filtered)
-
+        # Extract or compute velocity based on signal type
+        if self.signal_type == 'position':
+            # For position: compute velocity for SPARC
+            velocity = extract_velocity_from_position(signal, self.rate_hz)
+            # For jerk: use first dimension of position
+            signal_for_jerk = signal[:, 0] if signal.ndim > 1 else signal.squeeze()
+        else:  # self.signal_type == 'velocity'
+            # For velocity: use directly
+            if signal.ndim > 1 and signal.shape[1] > 1:
+                velocity = signal[:, 0]
+            else:
+                velocity = signal.squeeze()
+            signal_for_jerk = velocity
+
+        # Apply filtering to velocity for SPARC
+        filtered_velocity = self._filter_signal(velocity)
+        normalized = normalize_signal(filtered_velocity)
+
+        # Compute SPARC (always uses velocity)
         sparc = compute_sparc(normalized, self.rate_hz)
-        # Note: Since the smoothness module expects velocity input (as shown in examples), we specify signal_type='velocity' to compute proper jerk
-        jerk = compute_jerk_rms(filtered, self.rate_hz, signal_type='velocity')
+
+        # Compute jerk with appropriate signal type
+        filtered_for_jerk = self._filter_signal(signal_for_jerk)
+        jerk = compute_jerk_rms(filtered_for_jerk, self.rate_hz, signal_type=self.signal_type)
 
         return {"sparc": sparc, "jerk_rms": jerk}

pyeyesweb/utils/math_utils.py

@@ -203,4 +203,55 @@ def normalize_signal(signal):
         Returns original signal if max absolute value is 0.
     """
     max_val = np.max(np.abs(signal))
-    return signal / max_val if max_val != 0 else signal
+    return signal / max_val if max_val != 0 else signal
+
+
+def extract_velocity_from_position(position, rate_hz=50.0):
+    """Extract velocity from position data.
+
+    Computes velocity magnitude from position data of any dimensionality.
+    For 1D position, returns absolute velocity. For multi-dimensional position,
+    returns the Euclidean norm of the velocity vector.
+
+    Parameters
+    ----------
+    position : ndarray
+        Position data. Can be:
+        - 1D array: single position coordinate
+        - 2D array with shape (n_samples, n_dims): multi-dimensional positions
+    rate_hz : float, optional
+        Sampling rate in Hz (default: 50.0).
+
+    Returns
+    -------
+    ndarray
+        1D array of velocity magnitudes.
+
+    Examples
+    --------
+    >>> # 1D position data
+    >>> position_1d = np.array([0, 1, 2, 3, 4])
+    >>> velocity = extract_velocity_from_position(position_1d, rate_hz=100)
+
+    >>> # 2D position data (x, y coordinates)
+    >>> position_2d = np.array([[0, 0], [1, 0], [1, 1], [2, 1]])
+    >>> velocity = extract_velocity_from_position(position_2d, rate_hz=100)
+    """
+    rate_hz = validate_numeric(rate_hz, 'rate_hz', min_val=0.0001)
+    dt = 1.0 / rate_hz
+
+    position = np.asarray(position)
+
+    # Handle 1D position
+    if position.ndim == 1 or (position.ndim == 2 and position.shape[1] == 1):
+        signal_1d = position.squeeze()
+        return np.abs(np.gradient(signal_1d, dt))
+
+    # Handle multi-dimensional position
+    if position.ndim == 2:
+        # Compute derivatives along time axis (axis=0)
+        derivatives = np.gradient(position, dt, axis=0)
+        # Return Euclidean norm of velocity vector
+        return np.linalg.norm(derivatives, axis=1)
+
+    raise ValueError(f"Position must be 1D or 2D array, got shape {position.shape}")