ENH: Add framewise displacement peak detection function

Add framewise displacement peak detection function.

Take advantage of the commit to document the
`nifreeze.analysis.motion` module.

Co-authored-by: Oscar Esteban <[email protected]>

Author: jhlegarreta

src/nifreeze/analysis/filtering.py

@@ -0,0 +1,52 @@
+# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
+# vi: set ft=python sts=4 ts=4 sw=4 et:
+#
+# Copyright The NiPreps Developers <[email protected]>
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+# We support and encourage derived works from this project, please read
+# about our expectations at
+#
+#     https://www.nipreps.org/community/licensing/
+#
+"""Analysis data filtering."""
+
+import numpy as np
+
+
+def normalize(x: np.ndarray):
+    r"""Normalize data using the z-score.
+
+    The z-score normalization is computed as:
+
+    .. math::
+
+        z_i = \frac{x_i - \mu}{\sigma}
+
+    where $x_i$ is the framewise displacement at point $i$, $\mu$ is the mean
+    of all values, $\sigma$ is the standard deviation of the values, and $z_i$
+    is the normalized z-score.
+
+    Parameters
+    ----------
+    x : :obj:`~numpy.ndarray`
+        Data to be normalized.
+
+    Returns
+    -------
+    :obj:`~numpy.ndarray`
+        Normalized data.
+    """
+
+    return (x - np.mean(x)) / np.std(x)

src/nifreeze/analysis/motion.py

@@ -20,8 +20,10 @@
 #
 #     https://www.nipreps.org/community/licensing/
 #
+"""Motion analysis."""
 
 import numpy as np
+from scipy.stats import zscore
 
 
 def compute_percentage_change(
@@ -57,3 +59,33 @@ def compute_percentage_change(
     rel_diff[mask] = 100 * (test[mask] - reference[mask]) / reference[mask]
 
     return rel_diff
+
+
+def identify_spikes(fd: np.ndarray, threshold: float = 2.0):
+    """Identify motion spikes in framewise displacement data.
+
+    Identifies high-motion frames as timepoint exceeding a given threshold value
+    based on z-score normalized framewise displacement (FD) values.
+
+    Parameters
+    ----------
+    fd : :obj:`~numpy.ndarray`
+        Framewise displacement data.
+    threshold : :obj:`float`, optional
+        Threshold value to determine motion spikes.
+
+    Returns
+    -------
+    indices : :obj:`~numpy.ndarray`
+        Indices of identified motion spikes.
+    mask : :obj:`~numpy.ndarray`
+        Mask of identified motion spikes.
+    """
+
+    # Normalize (z-score)
+    fd_norm = zscore(fd)
+
+    mask = fd_norm > threshold
+    indices = np.where(mask)[0]
+
+    return indices, mask

test/test_analysis.py

@@ -31,6 +31,7 @@
     compute_z_score,
     identify_bland_altman_salient_data,
 )
+from nifreeze.analysis.motion import identify_spikes
 
 
 def test_compute_z_score():
@@ -141,3 +142,21 @@ def test_identify_bland_altman_salient_data():
 
     assert len(salient_data[BASalientEntity.RIGHT_INDICES.value]) == top_n
     assert len(salient_data[BASalientEntity.RIGHT_MASK.value]) == len(_data1)
+
+
+def test_identify_spikes(request):
+    rng = request.node.rng
+
+    n_samples = 450
+
+    fd = rng.normal(0, 5, n_samples)
+    threshold = 2.0
+
+    expected_indices = np.asarray([5, 57, 85, 100, 127, 180, 191, 202, 335, 393, 409])
+    expected_mask = np.zeros(n_samples, dtype=bool)
+    expected_mask[expected_indices] = True
+
+    obtained_indices, obtained_mask = identify_spikes(fd, threshold=threshold)
+
+    assert np.array_equal(obtained_indices, expected_indices)
+    assert np.array_equal(obtained_mask, expected_mask)