[FIX] Disable extrapolation of out-of-bounds volumes

CITATION.cff

@@ -674,3 +674,9 @@ authors:
     family-names: Nájera
     website: https://github.com/Titan-C
     affiliation: Checkmk
+  - given-names: Jordi
+    family-names: Huguet
+    website: https://github.com/jhuguetn
+    affiliation: BarcelonaBeta Brain Research Center
+    email: jhuguetn@gmail.com
+    orcid: https://orcid.org/0000-0001-8420-4833

doc/changes/latest.rst

@@ -23,3 +23,5 @@ Enhancements
 
 Changes
 -------
+
+- :bdg-success:`API` Expose scipy CubicSpline ``extrapolate`` parameter in :func:`~signal.clean` to control the interpolation of censored volumes in both ends of the BOLD signal data (:gh:`4028` by `Jordi Huguet`_).

nilearn/signal.py

@@ -547,6 +547,7 @@ def clean(
     high_pass=None,
     t_r=2.5,
     ensure_finite=False,
+    extrapolate=True,
     **kwargs,
 ):
     """Improve :term:`SNR` on masked :term:`fMRI` signals.
@@ -667,6 +668,11 @@ def clean(
         If `True`, the non-finite values (NANs and infs) found in the data
         will be replaced by zeros.
 
+    extrapolate : :obj:`bool`, default=True
+        If `True` and filter='butterworth', censored volumes in both ends of
+        the signal data will be interpolated before filtering. Otherwise, they
+        will be discarded from the band-pass filtering process.
+
     kwargs : dict
         Keyword arguments to be passed to functions called within ``clean``.
         Kwargs prefixed with ``'butterworth__'`` will be passed to
@@ -731,10 +737,9 @@ def clean(
         )
 
     # Interpolation / censoring
-    signals, confounds = _handle_scrubbed_volumes(
-        signals, confounds, sample_mask, filter_type, t_r
+    signals, confounds, sample_mask = _handle_scrubbed_volumes(
+        signals, confounds, sample_mask, filter_type, t_r, extrapolate
     )
-
     # Detrend
     # Detrend and filtering should apply to confounds, if confound presents
     # keep filters orthogonal (according to Lindquist et al. (2018))
@@ -810,19 +815,29 @@ def clean(
 
 
 def _handle_scrubbed_volumes(
-    signals, confounds, sample_mask, filter_type, t_r
+    signals, confounds, sample_mask, filter_type, t_r, extrapolate
 ):
     """Interpolate or censor scrubbed volumes."""
     if sample_mask is None:
-        return signals, confounds
+        return signals, confounds, sample_mask
 
     if filter_type == "butterworth":
-        signals = _interpolate_volumes(signals, sample_mask, t_r)
+        signals = _interpolate_volumes(signals, sample_mask, t_r, extrapolate)
+        # discard non-interpolated out-of-bounds volumes
+        signals = signals[~np.isnan(signals).all(axis=1), :]
         if confounds is not None:
-            confounds = _interpolate_volumes(confounds, sample_mask, t_r)
+            confounds = _interpolate_volumes(
+                confounds, sample_mask, t_r, extrapolate
+            )
+            # discard non-interpolated out-of-bounds volumes
+            confounds = confounds[~np.isnan(confounds).all(axis=1), :]
+        if sample_mask is not None and not extrapolate:
+            # reset the indexing of the sample_mask excluding non-interpolated
+            # volumes at the head of the data
+            sample_mask -= sample_mask[0]
     else:  # Or censor when no filtering, or cosine filter
         signals, confounds = _censor_signals(signals, confounds, sample_mask)
-    return signals, confounds
+    return signals, confounds, sample_mask
 
 
 def _censor_signals(signals, confounds, sample_mask):
@@ -833,12 +848,26 @@ def _censor_signals(signals, confounds, sample_mask):
     return signals, confounds
 
 
-def _interpolate_volumes(volumes, sample_mask, t_r):
+def _interpolate_volumes(volumes, sample_mask, t_r, extrapolate):
     """Interpolate censored volumes in signals/confounds."""
+    if extrapolate:
+        extrapolate_default = (
+            "By default the cubic spline interpolator extrapolates "
+            "the out-of-bounds censored volumes in the data run. This "
+            "can lead to undesired filtered signal results. Starting in "
+            "version 0.13, the default strategy will be not to extrapolate "
+            "but to discard those volumes at filtering."
+        )
+        warnings.warn(
+            category=FutureWarning,
+            message=extrapolate_default,
+        )
     frame_times = np.arange(volumes.shape[0]) * t_r
     remained_vol = frame_times[sample_mask]
     remained_x = volumes[sample_mask, :]
-    cubic_spline_fitter = CubicSpline(remained_vol, remained_x)
+    cubic_spline_fitter = CubicSpline(
+        remained_vol, remained_x, extrapolate=extrapolate
+    )
     volumes_interpolated = cubic_spline_fitter(frame_times)
     volumes[~sample_mask, :] = volumes_interpolated[~sample_mask, :]
     return volumes

nilearn/tests/test_signal.py

@@ -1065,8 +1065,9 @@ def test_handle_scrubbed_volumes():
     (
         interpolated_signals,
         interpolated_confounds,
+        sample_mask,
     ) = nisignal._handle_scrubbed_volumes(
-        signals, confounds, sample_mask, "butterworth", 2.5
+        signals, confounds, sample_mask, "butterworth", 2.5, True
     )
     np.testing.assert_equal(
         interpolated_signals[sample_mask, :], signals[sample_mask, :]
@@ -1075,8 +1076,79 @@ def test_handle_scrubbed_volumes():
         interpolated_confounds[sample_mask, :], confounds[sample_mask, :]
     )
 
-    scrubbed_signals, scrubbed_confounds = nisignal._handle_scrubbed_volumes(
-        signals, confounds, sample_mask, "cosine", 2.5
+    (
+        scrubbed_signals,
+        scrubbed_confounds,
+        sample_mask,
+    ) = nisignal._handle_scrubbed_volumes(
+        signals, confounds, sample_mask, "cosine", 2.5, True
     )
     np.testing.assert_equal(scrubbed_signals, signals[sample_mask, :])
     np.testing.assert_equal(scrubbed_confounds, confounds[sample_mask, :])
+
+
+def test_handle_scrubbed_volumes_extrapolation():
+    """Check interpolation of signals with and without extrapolation."""
+    signals, _, confounds = generate_signals(
+        n_features=11, n_confounds=5, length=40
+    )
+
+    censored_samples = 5
+    total_samples = censored_samples + 3
+    sample_mask = np.arange(signals.shape[0])
+    scrub_index = np.concatenate((np.arange(censored_samples), [10, 20, 30]))
+    sample_mask = np.delete(sample_mask, scrub_index)
+
+    # Test cubic spline interpolation (enabled extrapolation) in the
+    # very first n=5 samples of generated signal
+    extrapolate_warning = (
+        "By default the cubic spline interpolator extrapolates "
+        "the out-of-bounds censored volumes in the data run. This "
+        "can lead to undesired filtered signal results. Starting in "
+        "version 0.13, the default strategy will be not to extrapolate "
+        "but to discard those volumes at filtering."
+    )
+    with pytest.warns(FutureWarning, match=extrapolate_warning):
+        (
+            extrapolated_signals,
+            extrapolated_confounds,
+            extrapolated_sample_mask,
+        ) = nisignal._handle_scrubbed_volumes(
+            signals, confounds, sample_mask, "butterworth", 2.5, True
+        )
+    np.testing.assert_equal(signals.shape[0], extrapolated_signals.shape[0])
+    np.testing.assert_equal(
+        confounds.shape[0], extrapolated_confounds.shape[0]
+    )
+    np.testing.assert_equal(sample_mask, extrapolated_sample_mask)
+
+    # Test cubic spline interpolation without predicting values outside
+    # the range of the signal available (disabled extrapolation), discarding
+    # the first n censored samples of generated signal
+    (
+        interpolated_signals,
+        interpolated_confounds,
+        interpolated_sample_mask,
+    ) = nisignal._handle_scrubbed_volumes(
+        signals, confounds, sample_mask, "butterworth", 2.5, False
+    )
+    np.testing.assert_equal(
+        signals.shape[0], interpolated_signals.shape[0] + censored_samples
+    )
+    np.testing.assert_equal(
+        confounds.shape[0], interpolated_confounds.shape[0] + censored_samples
+    )
+    np.testing.assert_equal(sample_mask - sample_mask[0], interpolated_sample_mask)
+
+    # Assert that the modified sample mask (interpolated_sample_mask)
+    # can be applied to the interpolated signals and confounds
+    (
+        censored_signals,
+        censored_confounds,
+    ) = nisignal._censor_signals(interpolated_signals,
+                                 interpolated_confounds,
+                                 interpolated_sample_mask)
+    np.testing.assert_equal(
+        signals.shape[0], censored_signals.shape[0] + total_samples)
+    np.testing.assert_equal(
+        confounds.shape[0], censored_confounds.shape[0] + total_samples)