Add initial set of metrics

phys2denoise/due.py

@@ -0,0 +1,77 @@
+# emacs: at the end of the file
+# ex: set sts=4 ts=4 sw=4 et:
+# ## ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### #
+"""Stub file for a guaranteed safe import of duecredit constructs if duecredit is not available.
+
+To use it, place it into your project codebase to be imported, e.g. copy as
+
+    cp stub.py /path/tomodule/module/due.py
+
+Note that it might be better to avoid naming it duecredit.py to avoid shadowing
+installed duecredit.
+
+Then use in your code as
+
+    from .due import due, Doi, BibTeX, Text
+
+See  https://github.com/duecredit/duecredit/blob/master/README.md for examples.
+
+Origin:     Originally a part of the duecredit
+Copyright:  2015-2019  DueCredit developers
+License:    BSD-2
+"""
+
+__version__ = "0.0.8"
+
+
+class InactiveDueCreditCollector(object):
+    """Just a stub at the Collector which would not do anything."""
+
+    def _donothing(self, *args, **kwargs):
+        """Perform no good and no bad."""
+        pass
+
+    def dcite(self, *args, **kwargs):
+        """If I could cite I would."""
+
+        def nondecorating_decorator(func):
+            return func
+
+        return nondecorating_decorator
+
+    active = False
+    activate = add = cite = dump = load = _donothing
+
+    def __repr__(self):
+        """Perform magic."""
+        return self.__class__.__name__ + "()"
+
+
+def _donothing_func(*args, **kwargs):
+    """Perform no good and no bad."""
+    pass
+
+
+try:
+    from duecredit import BibTeX, Doi, Text, Url, due
+
+    if "due" in locals() and not hasattr(due, "cite"):
+        raise RuntimeError("Imported due lacks .cite. DueCredit is now disabled")
+except Exception as e:
+    if not isinstance(e, ImportError):
+        import logging
+
+        logging.getLogger("duecredit").error(
+            "Failed to import duecredit due to %s" % str(e)
+        )
+    # Initiate due stub
+    due = InactiveDueCreditCollector()
+    BibTeX = Doi = Url = Text = _donothing_func
+
+# Emacs mode definitions
+# Local Variables:
+# mode: python
+# py-indent-offset: 4
+# tab-width: 4
+# indent-tabs-mode: nil
+# End:

phys2denoise/metrics/__init__.py

@@ -0,0 +1,15 @@
+"""Metrics derived from physiological signals."""
+from .cardiac import crf, cardiac_phase
+from .chest_belt import rpv, env, rv, rrf, respiratory_phase
+from .multimodal import retroicor
+
+__all__ = [
+    "crf",
+    "cardiac_phase",
+    "rpv",
+    "env",
+    "rv",
+    "rrf",
+    "respiratory_phase",
+    "retroicor",
+]

phys2denoise/metrics/cardiac.py

@@ -0,0 +1,111 @@
+"""Denoising metrics for cardio recordings."""
+import numpy as np
+
+from .. import references
+from ..due import due
+
+
+def iht():
+    """Calculate instantaneous heart rate."""
+    pass
+
+
+@due.dcite(references.CHANG_GLOVER_2009)
+def crf(samplerate, oversampling=50, time_length=32, onset=0.0, tr=2.0):
+    """Calculate the cardiac response function using Chang and Glover's definition.
+
+    Parameters
+    ----------
+    samplerate : :obj:`float`
+        Sampling rate of data, in seconds.
+    oversampling : :obj:`int`, optional
+        Temporal oversampling factor, in seconds. Default is 50.
+    time_length : :obj:`int`, optional
+        RRF kernel length, in seconds. Default is 32.
+    onset : :obj:`float`, optional
+        Onset of the response, in seconds. Default is 0.
+
+    Returns
+    -------
+    crf : array-like
+        Cardiac or "heart" response function
+
+    Notes
+    -----
+    This cardiac response function was defined in [1]_, Appendix A.
+
+    The core code for this function comes from metco2, while several of the
+    parameters, including oversampling, time_length, and onset, are modeled on
+    nistats' HRF functions.
+
+    References
+    ----------
+    .. [1] C. Chang & G. H. Glover, "Relationship between respiration,
+       end-tidal CO2, and BOLD signals in resting-state fMRI," Neuroimage,
+       issue 47, vol. 4, pp. 1381-1393, 2009.
+    """
+
+    def _crf(t):
+        rf = 0.6 * t ** 2.7 * np.exp(-t / 1.6) - 16 * (
+            1 / np.sqrt(2 * np.pi * 9)
+        ) * np.exp(-0.5 * (((t - 12) ** 2) / 9))
+        return rf
+
+    dt = tr / oversampling
+    time_stamps = np.linspace(
+        0, time_length, np.rint(float(time_length) / dt).astype(np.int)
+    )
+    time_stamps -= onset
+    crf_arr = _crf(time_stamps)
+    crf_arr = crf_arr / max(abs(crf_arr))
+    return crf_arr
+
+
+def cardiac_phase(peaks, slice_timings, n_scans, t_r):
+    """Calculate cardiac phase from cardiac peaks.
+
+    Assumes that timing of cardiac events are given in same units
+    as slice timings, for example seconds.
+
+    Parameters
+    ----------
+    peaks : 1D array_like
+        Cardiac peak times, in seconds.
+    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.
+
+    Returns
+    -------
+    phase_card : array_like
+        Cardiac phase signal.
+    """
+    n_slices = np.shape(slice_timings)
+    phase_card = np.zeros(n_scans)
+
+    for i_slice in range(n_slices):
+        # find previous cardiac peaks
+        previous_card_peaks = np.asarray(np.nonzero(peaks < slice_timings[i_slice]))
+        if np.size(previous_card_peaks) == 0:
+            t1 = 0
+        else:
+            last_peak = previous_card_peaks[0][-1]
+            t1 = peaks[last_peak]
+
+        # find posterior cardiac peaks
+        next_card_peaks = np.asarray(np.nonzero(peaks > slice_timings[i_slice]))
+        if np.size(next_card_peaks) == 0:
+            t2 = n_scans * t_r
+        else:
+            next_peak = next_card_peaks[0][0]
+            t2 = peaks[next_peak]
+
+        # compute cardiac phase
+        phase_card[i_slice] = (2 * np.math.pi * (slice_timings[i_slice] - t1)) / (
+            t2 - t1
+        )
+
+    return phase_card