Complete overwrite of RETROICOR computation to allow multiple physiological modalities and their interactions

phys2denoise/metrics/__init__.py

@@ -1,9 +1,8 @@
 """Metrics derived from physiological signals."""
-from .cardiac import cardiac_phase, heart_beat_interval, heart_rate_variability
+from .cardiac import heart_beat_interval, heart_rate_variability
 from .chest_belt import (
     env,
     respiratory_pattern_variability,
-    respiratory_phase,
     respiratory_variance,
     respiratory_variance_time,
 )
@@ -12,8 +11,6 @@
 
 __all__ = [
     "retroicor",
-    "cardiac_phase",
-    "respiratory_phase",
     "heart_beat_interval",
     "heart_rate_variability",
     "env",

phys2denoise/metrics/cardiac.py

@@ -346,89 +346,3 @@ def heart_beat_interval(
     )
 
     return data, hbi
-
-
-@return_physio_or_metric()
-@physio.make_operation()
-def cardiac_phase(data, slice_timings, n_scans, t_r, peaks=None, fs=None, **kwargs):
-    """Calculate cardiac phase from cardiac peaks.
-
-    Assumes that timing of cardiac events are given in same units
-    as slice timings, for example seconds.
-
-    Parameters
-    ----------
-    data : physutils.Physio, np.ndarray, or array-like object
-        Object containing the timeseries of the recorded cardiac signal
-    slice_timings : 1D array_like
-        Slice times, in seconds.
-    n_scans : int
-        Number of volumes in the imaging run.
-    t_r : float
-        Sampling rate of the imaging run, in seconds.
-    fs : float, optional if data is a physutils.Physio
-        Sampling rate of `data` in Hz
-    peaks : :obj:`numpy.ndarray`, optional if data is a physutils.Physio
-        Indices of peaks in `data`
-
-    Returns
-    -------
-    phase_card : array_like
-        Cardiac phase signal, of shape (n_scans,)
-    """
-    if isinstance(data, physio.Physio):
-        # Initialize physio object
-        data = physio.check_physio(data, ensure_fs=True, copy=True)
-    elif fs is not None and peaks is not None:
-        data = io.load_physio(data, fs=fs)
-        data._metadata["peaks"] = peaks
-    else:
-        raise ValueError(
-            """
-            To use this function you should either provide a Physio object
-            with existing peaks metadata (e.g. using the peakdet module), or
-            by providing the physiological data timeseries, the sampling frequency,
-            and the peak indices separately.
-            """
-        )
-    if data.peaks.size == 0:
-        raise ValueError(
-            """
-            Peaks must be a non-empty list.
-            Make sure to run peak detection on your physiological data first,
-            using the peakdet module, or other software of your choice.
-            """
-        )
-
-    assert slice_timings.ndim == 1, "Slice times must be a 1D array"
-    n_slices = np.size(slice_timings)
-    phase_card = np.zeros((n_scans, n_slices))
-
-    card_peaks_sec = data.peaks / data.fs
-    for i_slice in range(n_slices):
-        # generate slice acquisition timings across all scans
-        times_crSlice = t_r * np.arange(n_scans) + slice_timings[i_slice]
-        phase_card_crSlice = np.zeros(n_scans)
-        for j_scan in range(n_scans):
-            previous_card_peaks = np.asarray(
-                np.nonzero(card_peaks_sec < times_crSlice[j_scan])
-            )
-            if np.size(previous_card_peaks) == 0:
-                t1 = 0
-            else:
-                last_peak = previous_card_peaks[0][-1]
-                t1 = card_peaks_sec[last_peak]
-            next_card_peaks = np.asarray(
-                np.nonzero(card_peaks_sec > times_crSlice[j_scan])
-            )
-            if np.size(next_card_peaks) == 0:
-                t2 = n_scans * t_r
-            else:
-                next_peak = next_card_peaks[0][0]
-                t2 = card_peaks_sec[next_peak]
-            phase_card_crSlice[j_scan] = (
-                2 * np.math.pi * (times_crSlice[j_scan] - t1)
-            ) / (t2 - t1)
-        phase_card[:, i_slice] = phase_card_crSlice
-
-    return data, phase_card

phys2denoise/metrics/chest_belt.py

@@ -306,72 +306,3 @@ def respiratory_variance(data, fs=None, window=6, **kwargs):
     rv_out = convolve_and_rescale(rv_arr, rrf(data.fs), rescale="zscore")
 
     return data, rv_out
-
-
-@return_physio_or_metric()
-@physio.make_operation()
-def respiratory_phase(data, n_scans, slice_timings, t_r, fs=None, **kwargs):
-    """Calculate respiratory phase from respiratory signal.
-
-    Parameters
-    ----------
-    data : physutils.Physio, np.ndarray, or array-like object
-        Object containing the timeseries of the recorded respiratory signal
-    n_scans
-        Number of volumes in the imaging run.
-    slice_timings
-        Slice times, in seconds.
-    t_r
-        Sample rate of the imaging run, in seconds.
-    fs : float, optional if data is a physutils.Physio
-        Sampling rate of `data` in Hz
-
-    Returns
-    -------
-    phase_resp : array_like
-        Respiratory phase signal, of shape (n_scans, n_slices).
-    """
-    if isinstance(data, physio.Physio):
-        # Initialize physio object
-        data = physio.check_physio(data, ensure_fs=True, copy=True)
-    elif fs is not None:
-        data = io.load_physio(data, fs=fs)
-    else:
-        raise ValueError(
-            """
-            To use this function you should either provide a Physio object
-            with the sampling frequency encapsulated, or
-            by providing the physiological data timeseries and the sampling
-            frequency separately.
-            """
-        )
-
-    assert slice_timings.ndim == 1, "Slice times must be a 1D array"
-    n_slices = np.size(slice_timings)
-    phase_resp = np.zeros((n_scans, n_slices))
-
-    # generate histogram from respiratory signal
-    # TODO: Replace with numpy.histogram
-    resp_hist, resp_hist_bins = np.histogram(data, bins=100)
-
-    # first compute derivative of respiration signal
-    resp_diff = np.diff(data, n=1)
-
-    for i_slice in range(n_slices):
-        # generate slice acquisition timings across all scans
-        times_crSlice = t_r * np.arange(n_scans) + slice_timings[i_slice]
-        phase_resp_crSlice = np.zeros(n_scans)
-        for j_scan in range(n_scans):
-            iphys = int(
-                max([1, round(times_crSlice[j_scan] * data.fs)])
-            )  # closest idx in resp waveform
-            iphys = min([iphys, len(resp_diff)])  # cannot be longer than resp_diff
-            thisBin = np.argmin(abs(data[iphys] - resp_hist_bins))
-            numerator = np.sum(resp_hist[0:thisBin])
-            phase_resp_crSlice[j_scan] = (
-                np.math.pi * np.sign(resp_diff[iphys]) * (numerator / len(data))
-            )
-
-        phase_resp[:, i_slice] = phase_resp_crSlice
-
-    return data, phase_resp