diff --git a/.all-contributorsrc b/.all-contributorsrc index de5f715..4e8b663 100644 --- a/.all-contributorsrc +++ b/.all-contributorsrc @@ -34,7 +34,8 @@ "infra", "projectManagement", "code", - "ideas" + "ideas", + "doc" ] }, { @@ -47,6 +48,15 @@ "ideas", "review" ] + }, + { + "login": "ineschh", + "name": "Inés Chavarría", + "avatar_url": "https://avatars.githubusercontent.com/u/72545702?v=4", + "profile": "https://github.com/ineschh", + "contributions": [ + "infra" + ] } ], "contributorsPerLine": 8, diff --git a/.github/workflows/pythonpublish.yml b/.github/workflows/pythonpublish.yml index 9e24640..4ddc1bd 100644 --- a/.github/workflows/pythonpublish.yml +++ b/.github/workflows/pythonpublish.yml @@ -18,7 +18,7 @@ jobs: - name: Set up Python uses: actions/setup-python@v1 with: - python-version: 3.7 + python-version: '3.7' - name: Install dependencies run: | python -m pip install --upgrade pip @@ -29,4 +29,4 @@ jobs: TWINE_PASSWORD: ${{ secrets.PYPI_PASSWORD }} run: | python -m python setup.py sdist bdist_wheel - python -m twine upload dist/* \ No newline at end of file + python -m twine upload dist/* diff --git a/CHANGELOG.md b/CHANGELOG.md new file mode 100644 index 0000000..126ffa6 --- /dev/null +++ b/CHANGELOG.md @@ -0,0 +1,47 @@ +# 0.2.0 (Thu Feb 25 2021) + +#### 💥 Breaking Change during development + +- Add initial set of metrics [#19](https://github.com/physiopy/phys2denoise/pull/19) ([@CesarCaballeroGaudes](https://github.com/CesarCaballeroGaudes) [@tsalo](https://github.com/tsalo) [@smoia](https://github.com/smoia)) + +#### Authors: 3 + +- Cesar Caballero Gaudes ([@CesarCaballeroGaudes](https://github.com/CesarCaballeroGaudes)) +- Stefano Moia ([@smoia](https://github.com/smoia)) +- Taylor Salo ([@tsalo](https://github.com/tsalo)) + +--- + +# 0.1.0 (Thu Feb 25 2021) + +:tada: This release contains work from new contributors! :tada: + +Thanks for all your work! + +:heart: Stefano Moia ([@smoia](https://github.com/smoia)) + +:heart: Inés Chavarría ([@ineschh](https://github.com/ineschh)) + +#### 💥 Breaking Change during development + +- Add main workflow of phys2denoise [#17](https://github.com/physiopy/phys2denoise/pull/17) ([@smoia](https://github.com/smoia) [@ineschh](https://github.com/ineschh)) + +#### ⚠️ Pushed to `master` + +- update all-contributors ([@smoia](https://github.com/smoia)) +- Correct attributes ([@smoia](https://github.com/smoia)) +- Add duecredit to setup ([@smoia](https://github.com/smoia)) +- General configuration of the infrastructure ([@smoia](https://github.com/smoia)) +- Correct development status ([@smoia](https://github.com/smoia)) +- Add infra files ([@smoia](https://github.com/smoia)) +- Initial commit ([@smoia](https://github.com/smoia)) + +#### 🏠 Internal + +- Update PR template to match `phys2bids` [#20](https://github.com/physiopy/phys2denoise/pull/20) ([@smoia](https://github.com/smoia)) +- Add initial infrastructure [#12](https://github.com/physiopy/phys2denoise/pull/12) ([@smoia](https://github.com/smoia)) + +#### Authors: 2 + +- Inés Chavarría ([@ineschh](https://github.com/ineschh)) +- Stefano Moia ([@smoia](https://github.com/smoia)) diff --git a/README.md b/README.md index 8742d16..d6b7b14 100644 --- a/README.md +++ b/README.md @@ -29,8 +29,9 @@ Thanks goes to these wonderful people ([emoji key](https://allcontributors.org/d
Cesar Caballero Gaudes
💻 🤔
Ross Markello
🚇 -
Stefano Moia
🔣 🚇 📆 💻 🤔 +
Stefano Moia
🔣 🚇 📆 💻 🤔 📖
Taylor Salo
💻 🤔 👀 +
Inés Chavarría
🚇 diff --git a/phys2denoise/cli/run.py b/phys2denoise/cli/run.py new file mode 100644 index 0000000..4103513 --- /dev/null +++ b/phys2denoise/cli/run.py @@ -0,0 +1,195 @@ +# -*- coding: utf-8 -*- +"""Parser for phys2denoise.""" + + +import argparse + +from phys2denoise import __version__ +from phys2denoise.metrics.cardiac import crf +from phys2denoise.metrics.chest_belt import rpv, rv, rvt, rrf, env + + +def _get_parser(): + """ + Parse command line inputs for this function. + + Returns + ------- + parser.parse_args() : argparse dict + + Notes + ----- + Default values must be updated in __call__ method from MetricsArgDict class. + # Argument parser follow template provided by RalphyZ. + # https://stackoverflow.com/a/43456577 + """ + + parser = argparse.ArgumentParser() + optional = parser._action_groups.pop() + required = parser.add_argument_group("Required Argument") + metrics = parser.add_argument_group("Metrics") + metric_arg = parser.add_argument_group("Metrics Arguments") + # Required arguments + required.add_argument("-in", "--input-file", + dest="filename", + type=str, + help="Full path and name of the file containing " + "physiological data, with or without extension.", + required=True) + # Important optional arguments + optional.add_argument("-outdir", "--output-dir", + dest="outdir", + type=str, + help="Folder where output should be placed. " + "Default is current folder.", + default=".") + # Metric selection + metrics.add_argument("-crf", "--cardiac-response-function", + dest="metrics", + action="append_const", + const=crf, + help="Cardiac response function. Requires the following " + "inputs:sample-rate, oversampling, time-length, " + "onset and tr.", + default=[]) + metrics.add_argument("-rpv", "--respiratory-pattern-variability", + dest="metrics", + action="append_const", + const=rpv, + help="Respiratory pattern variability. Requires the following " + "input: window.", + default=[]) + metrics.add_argument("-env", "--envelope", + dest="metrics", + action="append_const", + const=env, + help="Respiratory pattern variability calculated across a sliding " + "window. Requires the following inputs: sample-rate, window and lags.", + default=[]) + metrics.add_argument("-rv", "--respiratory-variance", + dest="metrics", + action="append_const", + const=rv, + help="Respiratory variance. Requires the following inputs: " + "sample-rate, window and lags. If the input file " + "not a .phys file, it also requires peaks and troughs", + default=[]) + """ + metrics.add_argument("-rvt", "--respiratory-volume-per-time", + dest="metrics", + action="append_const", + const="rvt", + help="Respiratory volume-per-time. Requires the following inputs: " + "sample-rate, window, lags, peaks and troughs.", + default=[]) + """ + metrics.add_argument("-rrf", "--respiratory-response-function", + dest="metrics", + action="append_const", + const=rrf, + help="Respiratory response function. Requires the following inputs: " + "sample-rate, oversampling, time-length, onset and tr.", + default=[]) + metrics.add_argument("-rcard", "--retroicor-card", + dest="metrics", + action="append_const", + const="r_card", + help="Computes regressors for cardiac signal. Requires the following " + "inputs: tr, nscans and n_harm.", + default=[]) + metrics.add_argument("-rresp", "--retroicor-resp", + dest="metrics", + action="append_const", + const="r_resp", + help="Computes regressors for respiratory signal. Requires the following " + "inputs: tr, nscans and n_harm.", + default=[]) + # Metric arguments + metric_arg.add_argument("-sr", "--sample-rate", + dest="sample_rate", + type=float, + help="Sampling rate of the physiological data in Hz.", + default=None) + metric_arg.add_argument("-pk", "--peaks", + dest="peaks", + type=str, + help="Full path and filename of the list with the indexed peaks' " + "positions of the physiological data.", + default=None) + metric_arg.add_argument("-tg", "--troughs", + dest="troughs", + type=str, + help="Full path and filename of the list with the indexed troughs' " + "positions of the physiological data.", + default=None) + metric_arg.add_argument("-os", "--oversampling", + dest="oversampling", + type=int, + help="Temporal oversampling factor. " + "Default is 50.", + default=50) + metric_arg.add_argument("-tl", "--time-length", + dest="time_length", + type=int, + help="RRF or CRF Kernel length in seconds.", + default=None) + metric_arg.add_argument("-onset", "--onset", + dest="onset", + type=float, + help="Onset of the response in seconds. " + "Default is 0.", + default=0) + metric_arg.add_argument("-tr", "--tr", + dest="tr", + type=float, + help="TR of sequence in seconds.", + default=None) + metric_arg.add_argument("-win", "--window", + dest="window", + type=int, + help="Size of the sliding window in seconds. " + "Default is 6 seconds.", + default=6) + metric_arg.add_argument("-lags", "--lags", + dest="lags", + nargs="*", + type=int, + help="List of lags to apply to the RV estimate " + "in seconds.", + default=None) + metric_arg.add_argument("-nscans", "--number-scans", + dest="nscans", + type=int, + help="Number of timepoints in the imaging data. " + "Also called sub-bricks, TRs, scans, volumes." + "Default is 1.", + default=1) + metric_arg.add_argument("-nharm", "--number-harmonics", + dest="n_harm", + type=int, + help="Number of harmonics.", + default=None) + + # Other optional arguments + optional.add_argument("-debug", "--debug", + dest="debug", + action="store_true", + help="Only print debugging info to log file. Default is False.", + default=False) + optional.add_argument("-quiet", "--quiet", + dest="quiet", + action="store_true", + help="Only print warnings to log file. Default is False.", + default=False) + optional.add_argument("-v", "--version", action="version", + version=("%(prog)s " + __version__)) + + parser._action_groups.append(optional) + + return parser + + +if __name__ == "__main__": + raise RuntimeError("phys2denoise/cli/run.py should not be run directly;\n" + "Please `pip install` phys2denoise and use the " + "`phys2denoise` command") diff --git a/phys2denoise/due.py b/phys2denoise/due.py new file mode 100644 index 0000000..4ff2457 --- /dev/null +++ b/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: diff --git a/phys2denoise/metrics/__init__.py b/phys2denoise/metrics/__init__.py index e69de29..8eb43c8 100644 --- a/phys2denoise/metrics/__init__.py +++ b/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", +] diff --git a/phys2denoise/metrics/cardiac.py b/phys2denoise/metrics/cardiac.py new file mode 100644 index 0000000..f52d7b2 --- /dev/null +++ b/phys2denoise/metrics/cardiac.py @@ -0,0 +1,119 @@ +"""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, sample_rate, 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. + sample_rate : float + Sample rate of physio, in Hertz. + 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, of shape (n_scans,) + """ + 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 = peaks / sample_rate + 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 phase_card diff --git a/phys2denoise/metrics/chest_belt.py b/phys2denoise/metrics/chest_belt.py index 1098bed..e6bb51e 100644 --- a/phys2denoise/metrics/chest_belt.py +++ b/phys2denoise/metrics/chest_belt.py @@ -1,7 +1,19 @@ """Denoising metrics for chest belt recordings.""" +import matplotlib as mpl import numpy as np +import pandas as pd +from scipy.ndimage.filters import convolve1d +from scipy.signal import detrend +from scipy.stats import zscore from scipy.interpolate import interp1d +mpl.use("TkAgg") +import matplotlib.pyplot as plt + +from .. import references +from ..due import due +from . import utils + def rvt(belt_ts, peaks, troughs, samplerate, lags=(0, 4, 8, 12)): """ @@ -63,3 +75,246 @@ def rvt(belt_ts, peaks, troughs, samplerate, lags=(0, 4, 8, 12)): rvt_lags[:, ind] = temp_rvt return rvt_lags + + +@due.dcite(references.POWER_2018) +def rpv(resp, window=6): + """Calculate respiratory pattern variability. + + Parameters + ---------- + resp : 1D array_like + window : int + + Returns + ------- + rpv_val : float + Respiratory pattern variability value. + + Notes + ----- + This metric was first introduced in [1]_. + + 1. Z-score respiratory belt signal + 2. Calculate upper envelope + 3. Calculate standard deviation of envelope + + References + ---------- + .. [1] J. D. Power et al., "Ridding fMRI data of motion-related influences: + Removal of signals with distinct spatial and physical bases in multiecho + data," Proceedings of the National Academy of Sciences, issue 9, vol. + 115, pp. 2105-2114, 2018. + """ + # First, z-score respiratory traces + resp_z = zscore(resp) + + # Collect upper envelope + rpv_upper_env = utils.rms_envelope_1d(resp_z, window) + + # Calculate standard deviation + rpv_val = np.std(rpv_upper_env) + return rpv_val + + +@due.dcite(references.POWER_2020) +def env(resp, samplerate, window=10): + """Calculate respiratory pattern variability across a sliding window. + + Parameters + ---------- + resp : (X,) :obj:`numpy.ndarray` + A 1D array with the respiratory belt time series. + samplerate : :obj:`float` + Sampling rate for resp, in Hertz. + window : :obj:`int`, optional + Size of the sliding window, in seconds. + Default is 10. + + Returns + ------- + env_arr + + Notes + ----- + This metric was first introduced in [1]_. + + Across a sliding window, do the following: + 1. Z-score respiratory belt signal + 2. Calculate upper envelope + 3. Calculate standard deviation of envelope + + References + ---------- + .. [1] J. D. Power et al., "Characteristics of respiratory measures in + young adults scanned at rest, including systematic changes and 'missed' + deep breaths," Neuroimage, vol. 204, 2020. + """ + # Convert window to Hertz + window = int(window * samplerate) + + # Calculate RPV across a rolling window + env_arr = ( + pd.Series(resp).rolling(window=window, center=True).apply(rpv, args=(window,)) + ) + env_arr[np.isnan(env_arr)] = 0.0 + return env_arr + + +@due.dcite(references.CHANG_GLOVER_2009) +def rv(resp, samplerate, window=6, lags=(0,)): + """Calculate respiratory variance. + + Parameters + ---------- + resp : (X,) :obj:`numpy.ndarray` + A 1D array with the respiratory belt time series. + samplerate : :obj:`float` + Sampling rate for resp, in Hertz. + window : :obj:`int`, optional + Size of the sliding window, in seconds. + Default is 6. + + Returns + ------- + rv_out : (X, 2) :obj:`numpy.ndarray` + Respiratory variance values. + The first column is raw RV values, after detrending/normalization. + The second column is RV values convolved with the RRF, after detrending/normalization. + + Notes + ----- + Respiratory variance (RV) was introduced in [1]_, and consists of the + standard deviation of the respiratory trace within a 6-second window. + + This metric is often lagged back and/or forward in time and convolved with + a respiratory response function before being included in a GLM. + Regressors also often have mean and linear trends removed and are + standardized prior to regressions. + + References + ---------- + .. [1] C. Chang & G. H. Glover, "Relationship between respiration, + end-tidal CO2, and BOLD signals in resting-state fMRI," Neuroimage, + issue 4, vol. 47, pp. 1381-1393, 2009. + """ + # Convert window to Hertz + window = int(window * samplerate) + + # Raw respiratory variance + rv_arr = pd.Series(resp).rolling(window=window, center=True).std() + rv_arr[np.isnan(rv_arr)] = 0.0 + + # Convolve with rrf + rrf_arr = rrf(samplerate, oversampling=1) + rv_convolved = convolve1d(rv_arr, rrf_arr, axis=0) + + # Concatenate the raw and convolved versions + rv_combined = np.stack((rv_arr, rv_convolved), axis=-1) + + # Detrend and normalize + rv_combined = rv_combined - np.mean(rv_combined, axis=0) + rv_combined = detrend(rv_combined, axis=0) + rv_out = zscore(rv_combined, axis=0) + return rv_out + + +@due.dcite(references.CHANG_GLOVER_2009) +def rrf(samplerate, oversampling=50, time_length=50, onset=0.0, tr=2.0): + """Calculate the respiratory 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. Default is 50. + time_length : :obj:`int`, optional + RRF kernel length, in seconds. Default is 50. + onset : :obj:`float`, optional + Onset of the response, in seconds. Default is 0. + + Returns + ------- + rrf : array-like + respiratory response function + + Notes + ----- + This respiratory 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 _rrf(t): + rf = 0.6 * t ** 2.1 * np.exp(-t / 1.6) - 0.0023 * t ** 3.54 * np.exp(-t / 4.25) + return rf + + dt = tr / oversampling + time_stamps = np.linspace( + 0, time_length, np.rint(float(time_length) / dt).astype(np.int) + ) + time_stamps -= onset + rrf_arr = _rrf(time_stamps) + rrf_arr = rrf_arr / max(abs(rrf_arr)) + return rrf_arr + + +def respiratory_phase(resp, sample_rate, n_scans, slice_timings, t_r): + """Calculate respiratory phase from respiratory signal. + + Parameters + ---------- + resp : 1D array_like + Respiratory signal. + sample_rate : float + Sample rate of physio, in Hertz. + 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. + + Returns + ------- + phase_resp : array_like + Respiratory phase signal, of shape (n_scans, n_slices). + """ + 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, _ = plt.hist(resp, bins=100) + + # first compute derivative of respiration signal + resp_diff = np.diff(resp, 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] * sample_rate)]) + ) # closest idx in resp waveform + iphys = min([iphys, len(resp_diff)]) # cannot be longer than resp_diff + thisBin = np.argmin(abs(resp[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(resp)) + ) + + phase_resp[:, i_slice] = phase_resp_crSlice + + return phase_resp \ No newline at end of file diff --git a/phys2denoise/metrics/gas.py b/phys2denoise/metrics/gas.py new file mode 100644 index 0000000..a4e3a59 --- /dev/null +++ b/phys2denoise/metrics/gas.py @@ -0,0 +1 @@ +"""Denoising metrics for gas recordings.""" diff --git a/phys2denoise/metrics/multimodal.py b/phys2denoise/metrics/multimodal.py new file mode 100644 index 0000000..fb95eae --- /dev/null +++ b/phys2denoise/metrics/multimodal.py @@ -0,0 +1,105 @@ +"""These functions compute RETROICOR regressors (Glover et al. 2000).""" + +import numpy as np + +from .cardiac import cardiac_phase +from .chest_belt import respiratory_phase +from .. import references +from ..due import due + + +@due.dcite(references.GLOVER_2000) +def retroicor( + physio, + sample_rate, + t_r, + n_scans, + slice_timings, + n_harmonics, + card=False, + resp=False, +): + """Compute RETROICOR regressors. + + Parameters + ---------- + physio : array_like + 1D array, whether cardiac or respiratory signal. + If cardiac, the array is a set of peaks in seconds. + If respiratory, the array is the actual respiratory signal. + sample_rate : float + Physio sample rate, in Hertz. + t_r : float + Imaging sample rate, in seconds. + n_scans : int + Number of volumes in the imaging run. + slice_timings : array_like + Slice times, in seconds. + n_harmonics : int + ??? + card : bool, optional + Whether the physio data correspond to cardiac or repiratory signal. + resp : bool, optional + Whether the physio data correspond to cardiac or repiratory signal. + + Returns + ------- + retroicor_regressors : list + phase : array_like + 2D array of shape (n_scans, n_slices) + + Notes + ----- + RETROICOR regressors should be regressed from the imaging data *before* + any other preprocessing, including slice-timing correction and motion correction. + + References + ---------- + * Glover, G. H., Li, T. Q., & Ress, D. (2000). + Image‐based method for retrospective correction of physiological + motion effects in fMRI: RETROICOR. + Magnetic Resonance in Medicine: + An Official Journal of the International Society for Magnetic Resonance in Medicine, + 44(1), 162-167. + """ + n_slices = np.shape(slice_timings) # number of slices + + # initialize output variables + retroicor_regressors = [] + phase = np.empty((n_scans, n_slices)) + + for i_slice in range(n_slices): + retroicor_regressors.append(np.empty((n_scans, 2 * n_harmonics))) + + # Initialize slice timings for current slice + crslice_timings = t_r * np.arange(n_scans) + slice_timings[i_slice] + + # Compute physiological phases using the timings of physio events (e.g. peaks) + # slice sampling times + if card: + phase[:, i_slice] = cardiac_phase( + physio, + crslice_timings, + n_scans, + t_r, + ) + + if resp: + phase[:, i_slice] = respiratory_phase( + physio, + sample_rate, + n_scans, + slice_timings, + t_r, + ) + + # Compute retroicor regressors + for j_harm in range(n_harmonics): + retroicor_regressors[i_slice][:, 2 * j_harm] = np.cos( + (j_harm + 1) * phase[i_slice] + ) + retroicor_regressors[i_slice][:, 2 * j_harm + 1] = np.sin( + (j_harm + 1) * phase[i_slice] + ) + + return retroicor_regressors, phase diff --git a/phys2denoise/metrics/utils.py b/phys2denoise/metrics/utils.py new file mode 100644 index 0000000..7426a5c --- /dev/null +++ b/phys2denoise/metrics/utils.py @@ -0,0 +1,125 @@ +"""Miscellaneous utility functions for metric calculation.""" +import logging + +import numpy as np + +LGR = logging.getLogger(__name__) +LGR.setLevel(logging.INFO) + + +def print_metric_call(metric, args): + """ + Log a message to describe how a metric is being called. + + Parameters + ---------- + metric : function + Metric function that is being called + args : dict + Dictionary containing all arguments that are used to parametrise metric + + Notes + ----- + Outcome + An info-level message for the logger. + """ + msg = f'The {metric} regressor will be computed using the following parameters:' + + for arg in args: + msg = f'{msg}\n {arg} = {args[arg]}' + + msg = f'{msg}\n' + + LGR.info(msg) + + +def mirrorpad_1d(arr, buffer=250): + """Pad both sides of array with flipped values from array of length 'buffer'. + + Parameters + ---------- + arr + buffer + + Returns + ------- + arr_out + """ + mirror = np.flip(arr, axis=0) + idx = range(arr.shape[0] - buffer, arr.shape[0]) + pre_mirror = np.take(mirror, idx, axis=0) + idx = range(0, buffer) + post_mirror = np.take(mirror, idx, axis=0) + arr_out = np.concatenate((pre_mirror, arr, post_mirror), axis=0) + return arr_out + + +def rms_envelope_1d(arr, window=500): + """Conceptual translation of MATLAB 2017b's envelope(X, x, 'rms') function. + + Parameters + ---------- + arr + window + + Returns + ------- + rms_env : numpy.ndarray + The upper envelope. + + Notes + ----- + https://www.mathworks.com/help/signal/ref/envelope.html + """ + assert arr.ndim == 1, "Input data must be 1D" + assert window % 2 == 0, "Window must be even" + n_t = arr.shape[0] + buf = int(window / 2) + + # Pad array at both ends + arr = np.copy(arr).astype(float) + mean = np.mean(arr) + arr -= mean + arr = mirrorpad_1d(arr, buffer=buf) + rms_env = np.empty(n_t) + for i in range(n_t): + # to match matlab + start_idx = i + buf + stop_idx = i + buf + window + + # but this is probably more appropriate + # start_idx = i + buf - 1 + # stop_idx = i + buf + window + window_arr = arr[start_idx:stop_idx] + rms = np.sqrt(np.mean(window_arr ** 2)) + rms_env[i] = rms + rms_env += mean + return rms_env + + +def apply_lags(arr1d, lags): + """Apply delays (lags) to an array. + + Parameters + ---------- + arr1d : (X,) :obj:`numpy.ndarray` + One-dimensional array to apply delays to. + lags : (Y,) :obj:`tuple` or :obj:`int` + Delays, in the same units as arr1d, to apply to arr1d. Can be negative, + zero, or positive integers. + + Returns + ------- + arr_with_lags : (X, Y) :obj:`numpy.ndarray` + arr1d shifted according to lags. Each column corresponds to a lag. + """ + arr_with_lags = np.zeros((arr1d.shape[0], len(lags))) + for i_lag, lag in enumerate(lags): + if lag < 0: + arr_delayed = np.hstack((arr1d[lag:], np.zeros(lag))) + elif lag > 0: + arr_delayed = np.hstack((np.zeros(lag), arr1d[lag:])) + else: + arr_delayed = arr1d.copy() + arr_with_lags[:, i_lag] = arr_delayed + return arr_with_lags diff --git a/phys2denoise/phys2denoise.py b/phys2denoise/phys2denoise.py new file mode 100644 index 0000000..0cc681e --- /dev/null +++ b/phys2denoise/phys2denoise.py @@ -0,0 +1,230 @@ +#!/usr/bin/env python3 + +""" +Phys2denoise is a python3 library meant to prepare physiological regressors for fMRI denoising. + +The project is under development. + +Copyright 2020, The physiopy community. +Please scroll to bottom to read full license. + +""" + +import datetime +import logging +import os +import sys +from inspect import signature, _empty + +import numpy as np +import pandas as pd + +from phys2denoise.cli.run import _get_parser +from phys2denoise.metrics.cardiac import crf +from phys2denoise.metrics.chest_belt import rpv, rv, rvt, rrf +from phys2denoise.metrics.retroicor import retroicor + +from . import __version__ +from .due import due, Doi + +LGR = logging.getLogger(__name__) +LGR.setLevel(logging.INFO) + + +def save_bash_call(outdir): + """ + Save the bash call into file `p2d_call.sh`. + + Parameters + ---------- + metric : function + Metric function to retrieve arguments for + metric_args : dict + Dictionary containing all arguments for all functions requested by the + user + """ + arg_str = ' '.join(sys.argv[1:]) + call_str = f'phys2denoise {arg_str}' + outdir = os.path.abspath(outdir) + log_path = os.path.join(outdir, 'code', 'logs') + os.makedirs(log_path) + f = open(os.path.join(log_path, 'p2d_call.sh'), "a") + f.write(f'#!bin/bash \n{call_str}') + f.close() + + +def select_input_args(metric, metric_args): + """ + Retrieve required args for metric from a dictionary of possible arguments. + + This function checks what parameters are accepted by a metric. + Then, for each parameter, check if the user provided it or not. + If they did not, but the parameter is required, throw an error - + unless it's "physio", reserved name for the timeseries input to a metric. + Otherwise, use the default. + + Parameters + ---------- + metric : function + Metric function to retrieve arguments for + metric_args : dict + Dictionary containing all arguments for all functions requested by the + user + + Returns + ------- + args : dict + Arguments to provide as input to metric + + Raises + ------ + ValueError + If a required argument is missing + + """ + args = {} + + # Check the parameters required by the metric and given by the user (see docstring) + for param in signature(metric).parameters.values(): + if param.name not in metric_args: + if param.default == _empty and param.name != 'physio': + raise ValueError(f'Missing parameter {param} required ' + f'to run {metric}') + else: + args[param.name] = param.default + else: + args[param.name] = metric_args[param.name] + + return args + + +@due.dcite( + Doi(''), + path='phys2denoise', + description='Creation of regressors for physiological denoising', + version=__version__, + cite_module=True) +def phys2denoise(filename, outdir='.', + metrics=[crf, rpv, rv, rvt, rrf, 'retroicor_card', 'retroicor_resp'], + debug=False, quiet=False, **kwargs): + """ + Run main workflow of phys2denoise. + + Runs the parser, does some checks on input, then computes the required metrics. + + Notes + ----- + Any metric argument should go into kwargs! + The code was greatly copied from phys2bids (copyright the physiopy community) + + """ + # Check options to make them internally coherent pt. I + # #!# This can probably be done while parsing? + outdir = os.path.abspath(outdir) + log_path = os.path.join(outdir, 'code', 'logs') + os.makedirs(log_path) + + # Create logfile name + basename = 'phys2denoise_' + extension = 'tsv' + isotime = datetime.datetime.now().strftime('%Y-%m-%dT%H%M%S') + logname = os.path.join(log_path, (basename + isotime + '.' + extension)) + + # Set logging format + log_formatter = logging.Formatter( + '%(asctime)s\t%(name)-12s\t%(levelname)-8s\t%(message)s', + datefmt='%Y-%m-%dT%H:%M:%S') + + # Set up logging file and open it for writing + log_handler = logging.FileHandler(logname) + log_handler.setFormatter(log_formatter) + sh = logging.StreamHandler() + + if quiet: + logging.basicConfig(level=logging.WARNING, + handlers=[log_handler, sh], + format='%(levelname)-10s %(message)s') + elif debug: + logging.basicConfig(level=logging.DEBUG, + handlers=[log_handler, sh], + format='%(levelname)-10s %(message)s') + else: + logging.basicConfig(level=logging.INFO, + handlers=[log_handler, sh], + format='%(levelname)-10s %(message)s') + + version_number = __version__ + LGR.info(f'Currently running phys2denoise version {version_number}') + LGR.info(f'Input file is {filename}') + + # Check options to make them internally coherent pt. II + # #!# This can probably be done while parsing? + # filename, ftype = utils.check_input_type(filename) + + if not os.path.isfile(filename) and filename is not None: + raise FileNotFoundError(f'The file {filename} does not exist!') + + # Read input file + physio = np.genfromtxt(filename) + + # Prepare pandas dataset + regr = pd.DataFrame() + + # Goes through the list of metrics and calls them + for metric in metrics: + if metric == 'retroicor_card': + args = select_input_args(retroicor, kwargs) + args['card'] = True + retroicor_regrs = retroicor(physio, **args) + for vslice in range(len(args['slice_timings'])): + for harm in range(args['n_harm']): + key = f'rcor-card_s-{vslice}_hrm-{harm}' + regr[f'{key}_cos'] = retroicor_regrs[vslice][:, harm*2] + regr[f'{key}_sin'] = retroicor_regrs[vslice][:, harm*2+1] + elif metric == 'retroicor_resp': + args = select_input_args(retroicor, kwargs) + args['resp'] = True + retroicor_regrs = retroicor(physio, **args) + for vslice in range(len(args['slice_timings'])): + for harm in range(args['n_harm']): + key = f'rcor-resp_s-{vslice}_hrm-{harm}' + regr[f'{key}_cos'] = retroicor_regrs[vslice][:, harm*2] + regr[f'{key}_sin'] = retroicor_regrs[vslice][:, harm*2+1] + else: + args = select_input_args(metric, kwargs) + regr[metric.__name__] = metric(physio, **args) + + #!# Add regressors visualisation + + # Export regressors and sidecar + out_filename = os.join(outdir, 'derivatives', filename) + regr.to_csv(out_filename, sep='\t', index=False, float_format='%.6e') + #!# Add sidecar export + + +def _main(argv=None): + options = _get_parser().parse_args(argv) + + save_bash_call(options.outdir) + + phys2denoise(**vars(options)) + + +if __name__ == '__main__': + _main(sys.argv[1:]) + +""" +Copyright 2019, The phys2denoise community. + +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. +""" diff --git a/phys2denoise/references.py b/phys2denoise/references.py new file mode 100644 index 0000000..eeb914a --- /dev/null +++ b/phys2denoise/references.py @@ -0,0 +1,10 @@ +"""References to be imported and injected throughout the package.""" +from .due import Doi + +CHANG_GLOVER_2009 = Doi("10.1016/j.neuroimage.2009.04.048") + +POWER_2018 = Doi("10.1073/pnas.1720985115") + +POWER_2020 = Doi("10.1016/j.neuroimage.2019.116234") + +GLOVER_2000 = Doi("10.1002/1522-2594(200007)44:1<162::AID-MRM23>3.0.CO;2-E") diff --git a/phys2denoise/tests/test_metrics_cardiac.py b/phys2denoise/tests/test_metrics_cardiac.py new file mode 100644 index 0000000..bb10d76 --- /dev/null +++ b/phys2denoise/tests/test_metrics_cardiac.py @@ -0,0 +1,41 @@ +"""Tests for phys2denoise.metrics.cardiac.""" +import numpy as np + +from phys2denoise.metrics import cardiac + + +def test_crf_smoke(): + """Basic smoke test for CRF calculation.""" + samplerate = 0.01 # in seconds + oversampling = 20 + time_length = 20 + onset = 0 + tr = 0.72 + crf_arr = cardiac.crf( + samplerate, + oversampling=oversampling, + time_length=time_length, + onset=onset, + tr=tr, + ) + pred_len = np.rint(time_length / (tr / oversampling)) + assert crf_arr.ndim == 1 + assert crf_arr.shape == pred_len + + +def test_cardiac_phase_smoke(): + """Basic smoke test for cardiac phase calculation.""" + t_r = 1.0 + n_scans = 200 + sample_rate = 1 / 0.01 + slice_timings = np.linspace(0, t_r, 22)[1:-1] + peaks = np.array([0.534, 0.577, 10.45, 20.66, 50.55, 90.22]) + card_phase = cardiac.cardiac_phase( + peaks, + sample_rate=sample_rate, + slice_timings=slice_timings, + n_scans=n_scans, + t_r=t_r, + ) + assert card_phase.ndim == 2 + assert card_phase.shape == (n_scans, slice_timings.size) diff --git a/phys2denoise/tests/test_metrics_chest_belt.py b/phys2denoise/tests/test_metrics_chest_belt.py new file mode 100644 index 0000000..f256918 --- /dev/null +++ b/phys2denoise/tests/test_metrics_chest_belt.py @@ -0,0 +1,73 @@ +"""Tests for phys2denoise.metrics.chest_belt.""" +import numpy as np + +from phys2denoise.metrics import chest_belt + + +def test_rrf_smoke(): + """Basic smoke test for RRF calculation.""" + samplerate = 0.01 # in seconds + oversampling = 20 + time_length = 20 + onset = 0 + tr = 0.72 + rrf_arr = chest_belt.rrf( + samplerate, + oversampling=oversampling, + time_length=time_length, + onset=onset, + tr=tr, + ) + pred_len = int(np.rint(time_length / (tr / oversampling))) + assert rrf_arr.ndim == 1 + assert rrf_arr.size == pred_len + + +def test_respiratory_phase_smoke(): + """Basic smoke test for respiratory phase calculation.""" + t_r = 1.0 + n_scans = 200 + sample_rate = 1 / 0.01 + slice_timings = np.linspace(0, t_r, 22)[1:-1] + n_samples = int(np.rint((n_scans * t_r) * sample_rate)) + resp = np.random.normal(size=n_samples) + resp_phase = chest_belt.respiratory_phase( + resp, + sample_rate=sample_rate, + slice_timings=slice_timings, + n_scans=n_scans, + t_r=t_r, + ) + assert resp_phase.ndim == 2 + assert resp_phase.shape == (n_scans, slice_timings.size) + + +def test_rpv_smoke(): + """Basic smoke test for respiratory pattern variability calculation.""" + n_samples = 2000 + resp = np.random.normal(size=n_samples) + window = 50 + rpv_val = chest_belt.rpv(resp, window) + assert isinstance(rpv_val, float) + + +def test_env_smoke(): + """Basic smoke test for ENV calculation.""" + n_samples = 2000 + resp = np.random.normal(size=n_samples) + samplerate = 1 / 0.01 + window = 6 + env_arr = chest_belt.env(resp, samplerate=samplerate, window=window) + assert env_arr.ndim == 1 + assert env_arr.shape == (n_samples,) + + +def test_rv_smoke(): + """Basic smoke test for respiratory variance calculation.""" + n_samples = 2000 + resp = np.random.normal(size=n_samples) + samplerate = 1 / 0.01 + window = 6 + rv_arr = chest_belt.rv(resp, samplerate=samplerate, window=window) + assert rv_arr.ndim == 2 + assert rv_arr.shape == (n_samples, 2)