Merge branch 'dev'

* dev:
  interface complete
  Update retroicor.py
  Update README.md
  draft some utility functions
  draft the CLI

Author: htwangtw

README.md

@@ -1,11 +1,20 @@
 # PyRETROICOR
 Python implementation of RETROICOR
 
-## Dependency
-Python 3  
-numpy >= 1.18
+## Dependencies
+**[Python 3](https://www.python.org/)**
+| Package                                                    | Tested version |
+|------------------------------------------------------------|----------------|
+| [NumPy](http://www.numpy.org/)                             | 1.18.4         |
+| [Nibabel](https://nipy.org/nibabel/)                       | 3.1.0          |
 
 ## Reference
 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.
 
 MATLAB implementation from: https://github.com/mkassinopoulos/PRF_estimation
+
+## Support
+Please use [GitHub issues](https://github.com/htwangtw/pyretroicor/issues) for questions, bug reports or feature requests.
+
+## License
+This project is licensed under [BSD-3-Clause](https://opensource.org/licenses/BSD-3-Clause).

pyretroicor/__init__.py

@@ -0,0 +1,4 @@
+import pkg_resources
+
+__version__ = pkg_resources.require(
+    "pyretroicor")[0].version

pyretroicor/__main__.py

@@ -0,0 +1,109 @@
+"""
+Main entry point
+"""
+from pathlib import Path
+import argparse
+import json
+import gzip
+
+import numpy as np
+import nibabel as nb
+
+from pyretroicor.retroicor import retroicor_cardiac, retroicor_respiratory
+from pyretroicor.prepro import process_resp, process_cardiac
+
+def main():
+    parser = argparse.ArgumentParser()
+
+    parser.add_argument(
+        'filename',  metavar='path',
+        help="Path to fMRI nifti files."
+        )
+    parser.add_argument(
+        'physio',  metavar='physio',
+        help="Path to physiology files in BIDS standard. \
+        Accept cardiac and/or respiratory recording."
+        )
+    parser.add_argument(
+        'output',  metavar='output',
+        help="Path to save calculated derivatives."
+        )
+
+    args = parser.parse_args()
+    NII_FILE = args.filename
+    PHYSIO_FILE = args.physio
+    OUTDIR = args.output
+
+    # check all files neesed exist: nifti, physio, accompanied jsons
+    physio_json = PHYSIO_FILE.split(".tsv.gz")[0] + ".json"
+    func_json = NII_FILE.split(".nii.gz")[0] + ".json"
+    with open(func_json) as f:
+        func_meta = json.load(f) 
+
+    if not Path(physio_json).is_file():
+        raise ValueError("Missing json file for physiology data. Is data in BIDS format?")
+    else:
+        with open(physio_json) as f:
+            physio_meta = json.load(f) 
+
+        sampling_freq = physio_meta["SamplingFrequency"]
+        start_time = physio_meta["StartTime"]
+        columns = physio_meta["Columns"]
+
+
+    func = nb.load(NII_FILE)  # nibabel should do the test itself
+    physio = np.loadtxt(PHYSIO_FILE) 
+    # create time
+    time = np.arange(0, physio.shape[0])  * (1 / sampling_freq) + start_time
+    if start_time < 0:
+        time_zero = np.where(time < 0)[0][-1] + 1
+        time = time[time_zero:]
+        physio = physio[time_zero:, :]
+
+    cardiac = physio[:, columns.index("cardiac")]
+    respiratory = physio[:, columns.index("respiratory")]
+    
+    # TR labels
+    n_vol = func.shape[0]
+    tr = func_meta["RepetitionTime"]
+    tr_time = np.arange(1, n_vol + 1)  * tr
+    tr_idx = []
+    for t in tr_time:
+        idx = np.where(time < t)[0][-1] + 1
+        tr_idx.append(idx)
+
+    # peak detection and cleaning
+    r_peak = process_cardiac(cardiac, sampling_freq)  # r peak index
+    r_peak = time[r_peak]  # r peak time point
+    retro_card = retroicor_cardiac(time, r_peak)
+
+    # run retroicor
+    respiratory = process_resp(respiratory, sampling_freq)
+    retro_resp = retroicor_respiratory(respiratory)
+    # rextract value per TR
+    retroicor_regressors = np.hstack((retro_card, retro_resp))[tr_idx, :]
+
+    # create output dir
+    nii_file = Path(NII_FILE)
+    subject = nii_file.name.split("_")[0].split("sub-")[-1]
+    filename_root = nii_file.name.split("_bold")[0]
+    if "ses" in NII_FILE:
+        session = nii_file.name.split("_")[1].split("ses-")[-1]
+        filepath = Path(OUTDIR) / "pyretroicor" / f"sub-{subject}" / f"ses-{session}" 
+        out_name = f"{filename_root}_desc-retroicor_regressors.tsv"
+    else:
+        filepath = Path(OUTDIR) / "pyretroicor" / f"sub-{subject}"
+        out_name = f"{filename_root}_desc-retroicor_regressors.tsv"
+
+    filepath.mkdir(parents=True, exist_ok=True)
+
+    # save results as tsv
+    header = [f"retroicor_cardiac_{i + 1}" for i in range(6)] \
+        + [f"retroicor_respiration_{i + 1}" for i in range(6)]
+
+    np.savetxt(filepath / out_name, retroicor_regressors, 
+               delimiter="\t", header="\t".join(header), 
+               comments="")
+
+if __name__ == "__main__":
+    main()

pyretroicor/prepro.py

@@ -0,0 +1,101 @@
+"""physio data preprocessing"""
+
+import numpy as np
+
+def process_cardiac(cardiac, fs):
+    '''
+    preprocess cardiac signal and get R peaks times
+    '''
+    # clean outliers
+    cardiac = smooth(cardiac, fs)
+    outliers = rolling_mad(cardiac, int(0.5 * fs))
+    time = np.arange(0, cardiac.shape[0])
+    clean = np.interp(time, time[~outliers], cardiac[~outliers])
+    # get peaks
+    peak_idx = detect_peaks(clean, mpd=int(0.5 * fs))
+    return peak_idx
+
+def zscore(x):
+    m = np.mean(x)
+    sd = np.std(x)
+    return (x - m) / sd
+
+def detect_peaks(x, mpd):
+    # standardlize signal 
+    x = zscore(x)
+    dx = np.diff(x)
+    # find the edges of raises (proximate peaks)
+    idx = np.where((np.hstack((dx, 0)) <= 0) & (np.hstack((0, dx)) > 0))[0]
+    if idx[0] == 0:  # first and least location must not be peak
+        idx = idx[1:]
+    if idx[-1] == (x.size - 1):
+        idx = idx[:-1]
+
+    # remove peak smaller than average signal
+    idx = idx[x[idx] >= 0]
+
+    # remove peaks closer than mpd
+    idx = idx[np.argsort(x[idx])][::-1]  # sort idx by peak height
+    idel = np.zeros(idx.size, dtype=bool)
+    for i in range(idx.size):
+        if not idel[i]:
+            idel = idel | (idx >= idx[i] - mpd) & (idx <= idx[i] + mpd)
+            idel[i] = 0  # Keep current peak
+    # remove the small peaks and sort back the indices by their occurrence
+    idx = np.sort(idx[~idel])
+    return idx
+
+
+def process_resp(resp, fs):
+    '''
+    process respiratory data
+    '''
+    resp_f = smooth(resp, fs)
+    resp_der = np.diff(resp_f)
+    resp_der = np.append(resp_der, resp_der[-1])
+    time = np.arange(0, resp_der.shape[0])
+    outliers = rolling_mad(resp_der, int(0.5 * fs))
+    clean = np.interp(time, time[~outliers], resp_der[~outliers])
+    return clean
+
+def mad(arr):
+    """ Median Absolute Deviation: a "Robust" version of standard deviation.
+        Indices variabililty of the sample.
+        https://en.wikipedia.org/wiki/Median_absolute_deviation
+    """
+    med = np.median(arr)
+    return med, np.median(np.abs(arr - med))
+
+def rolling_mad(a, win):
+    '''
+    Rolling MAD outlier detection
+    '''
+    outliers = []
+    for i in range(win, len(a)):
+        cur = a[(i - win):i]
+        med, cur_mad = mad(cur)
+        cur_out = cur > (med + cur_mad * 3)
+        idx = list(np.arange((i - win), i)[cur_out])
+        outliers += idx
+    outliers = list(set(outliers))
+
+    # turn index into boolean
+    bool_outliers = np.zeros(a.shape[0], dtype=bool)
+    bool_outliers[outliers] = True
+    return bool_outliers
+
+
+def smooth(a, winsize):
+    '''
+    a: NumPy 1-D array containing the data to be smoothed
+    winsize: smoothing window size needs, which must be odd number,
+            as in the original MATLAB implementation
+    '''
+    winsize = int(winsize)
+    if winsize % 2 != 1:
+        winsize += 1
+    out0 = np.convolve(a, np.ones(winsize, dtype=int),'valid') / winsize
+    r = np.arange(1, winsize - 1, 2)
+    start = np.cumsum(a[: winsize - 1])[::2] / r
+    stop = (np.cumsum(a[: -winsize:-1])[::2] / r)[::-1]
+    return np.concatenate((start, out0, stop))

pyretroicor/retroicor.py

@@ -0,0 +1,71 @@
+import numpy as np
+
+def retroicor_cardiac(time, r_peak, M=3):
+    '''
+    we need documentations
+    translated from matlab code
+
+    r_peak - time of the r peaks
+    time - tr time stamp?
+    M - ??????
+    '''
+
+    # the algorithm
+    num_time = len(time)
+    cardiac_phase = np.zeros(num_time)
+    for i, t in enumerate(time):
+        closest_peak = np.abs(r_peak - t)
+        peak_idx = np.where(closest_peak == np.amin(closest_peak))[0][0]
+        min_left = (t - r_peak[peak_idx]) > 0  # check if the peak is the first or last
+
+        if peak_idx == 0 and not min_left:
+            t2 = r_peak[peak_idx]
+            t1 = t2 - 1
+        elif peak_idx == (len(r_peak) - 1) and min_left:
+            t1 = r_peak[peak_idx]
+            t2 = t1 + 1
+        elif min_left:
+            t1 = r_peak[peak_idx]
+            t2 = r_peak[peak_idx + 1]
+        else:
+            t1 = r_peak[peak_idx - 1]
+            t2 = r_peak[peak_idx]
+        cardiac_phase[i] = 2 * np.pi * (t - t1) / (t2 - t1) # Eq 2
+
+    regr = np.zeros((num_time, M * 2))
+    for i in range(M):
+        regr[:,i * 2] = np.cos(i * cardiac_phase)
+        regr[:,(i * 2) + 1] = np.sin(i * cardiac_phase)
+
+    return regr
+
+
+def retroicor_respiratory(resp, M=3):
+    '''
+    need doc
+    '''
+    num_time = len(resp)
+    resp_phase = np.zeros(num_time)
+
+    # the algorithm
+    nb = 500
+    val, edges = np.histogram(resp, bins=nb)
+
+    for i in range(num_time):
+        v = resp[i]
+        closest_peak = np.abs(v - edges)
+        edge = np.where(closest_peak == np.amin(closest_peak))[0][0]
+        if (edge - v) > 0:
+            edge = edge - 1
+        area = np.sum(val[0: edge])
+        sign_resp = np.sign(resp[i])
+        if sign_resp == 0:
+            sign_resp = 1
+        resp_phase[i] = np.pi * area * sign_resp / num_time
+
+    regr = np.zeros((num_time, M * 2))
+    for i in range(M):
+        regr[:,i * 2] = np.cos(i * resp_phase)
+        regr[:,(i * 2) + 1] = np.sin(i * resp_phase)
+
+    return regr

pyretroicor/util.py

@@ -0,0 +1 @@
+"""I/O and file path senity check"""

requirements.txt

@@ -0,0 +1,2 @@
+nibabel>=2.5
+numpy>=1.18