diff --git a/README.md b/README.md index 6ac70b93f..2872b70e8 100644 --- a/README.md +++ b/README.md @@ -21,7 +21,7 @@ For a review, see [Kundu et al. (2017), _NeuroImage_](https://paperpile.com/shar In tedana, we combine all collected echos, then decompose the resulting time series into components that can be classified as BOLD or non-BOLD based. This is performed in a series of steps including: -* Principle components analysis +* Principal components analysis * Independent components analysis * Component classification diff --git a/tedana/interfaces/t2smap.py b/tedana/interfaces/t2smap.py index bfbe67e76..067715978 100644 --- a/tedana/interfaces/t2smap.py +++ b/tedana/interfaces/t2smap.py @@ -1,54 +1,124 @@ +""" +""" import numpy as np import nibabel as nib -from tedana.utils import (niwrite, cat2echos, - makeadmask, unmask, fmask) +from tedana.utils import (niwrite, cat2echos, makeadmask, unmask, fmask) -def t2sadmap(catd, mask, tes, masksum, start_echo): +def fit(data, mask, tes, masksum, start_echo): """ - t2sadmap(catd,mask,tes,masksum) - - Input: + Fit voxel- and timepoint-wise monoexponential decay models to estimate + T2* and S0 timeseries. + """ + nx, ny, nz, n_echoes, n_trs = data.shape + echodata = fmask(data, mask) + tes = np.array(tes) - catd has shape (nx,ny,nz,Ne,nt) - mask has shape (nx,ny,nz) - tes is a 1d numpy array - masksum + t2sa_ts = np.zeros([nx, ny, nz, n_trs]) + s0va_ts = np.zeros([nx, ny, nz, n_trs]) + t2saf_ts = np.zeros([nx, ny, nz, n_trs]) + s0vaf_ts = np.zeros([nx, ny, nz, n_trs]) + + for vol in range(echodata.shape[-1]): + t2ss = np.zeros([nx, ny, nz, n_echoes - 1]) + s0vs = t2ss.copy() + # Fit monoexponential decay first for first echo only, + # then first two echoes, etc. + for i_echo in range(start_echo, n_echoes + 1): + B = np.abs(echodata[:, :i_echo, vol]) + 1 + B = np.log(B).transpose() + neg_tes = -1 * tes[:i_echo] + + # First row is constant, second is TEs for decay curve + # Independent variables for least-squares model + x = np.array([np.ones(i_echo), neg_tes]) + X = np.sort(x)[:, ::-1].transpose() + + beta, _, _, _ = np.linalg.lstsq(X, B) + t2s = 1. / beta[1, :].transpose() + s0 = np.exp(beta[0, :]).transpose() + + t2s[np.isinf(t2s)] = 500. + s0[np.isnan(s0)] = 0. + + t2ss[:, :, :, i_echo-2] = np.squeeze(unmask(t2s, mask)) + s0vs[:, :, :, i_echo-2] = np.squeeze(unmask(s0, mask)) + + # Limited T2* and S0 maps + fl = np.zeros([nx, ny, nz, len(tes)-1], bool) + for i_echo in range(n_echoes - 1): + fl_ = np.squeeze(fl[:, :, :, i_echo]) + fl_[masksum == i_echo + 2] = True + fl[:, :, :, i_echo] = fl_ + t2sa = np.squeeze(unmask(t2ss[fl], masksum > 1)) + s0va = np.squeeze(unmask(s0vs[fl], masksum > 1)) + + # Full T2* maps with S0 estimation errors + t2saf = t2sa.copy() + s0vaf = s0va.copy() + t2saf[masksum == 1] = t2ss[masksum == 1, 0] + s0vaf[masksum == 1] = s0vs[masksum == 1, 0] + + t2sa_ts[:, :, :, vol] = t2sa + s0va_ts[:, :, :, vol] = s0va + t2saf_ts[:, :, :, vol] = t2saf + s0vaf_ts[:, :, :, vol] = s0vaf + + return t2sa_ts, s0va_ts, t2saf_ts, s0vaf_ts + + +def t2sadmap(data, mask, tes, masksum, start_echo): """ - nx, ny, nz, Ne, nt = catd.shape - echodata = fmask(catd, mask) - Nm = echodata.shape[0] + Fit voxelwise monoexponential decay models to estimate T2* and S0 maps. + t2sadmap(data,mask,tes,masksum) - t2ss = np.zeros([nx, ny, nz, Ne - 1]) + Input: + data : :obj:`numpy.ndarray` + Concatenated BOLD data. Has shape (nx, ny, nz, n_echoes, n_trs) + mask : :obj:`numpy.ndarray` + Brain mask in 3D array. Has shape (nx, ny, nz) + tes : :obj:`numpy.ndarray` + Array of TEs, in milliseconds. + masksum : + """ + nx, ny, nz, n_echoes, n_trs = data.shape + echodata = fmask(data, mask) + n_voxels = echodata.shape[0] + tes = np.array(tes) + t2ss = np.zeros([nx, ny, nz, n_echoes - 1]) s0vs = t2ss.copy() - for ne in range(start_echo, Ne + 1): - + # Fit monoexponential decay first for first echo only, + # then first two echoes, etc. + for i_echo in range(start_echo, n_echoes + 1): # Do Log Linear fit - B = np.reshape(np.abs(echodata[:, :ne]) + 1, (Nm, ne * nt)).transpose() + B = np.reshape(np.abs(echodata[:, :i_echo, :]) + 1, + (n_voxels, i_echo*n_trs)).transpose() B = np.log(B) - neg_tes = [-1 * te for te in tes[:ne]] - x = np.array([np.ones(ne), neg_tes]) - X = np.tile(x, (1, nt)) + neg_tes = -1 * tes[:i_echo] + + # First row is constant, second is TEs for decay curve + # Independent variables for least-squares model + x = np.array([np.ones(i_echo), neg_tes]) + X = np.tile(x, (1, n_trs)) X = np.sort(X)[:, ::-1].transpose() - beta, res, rank, sing = np.linalg.lstsq(X, B) - t2s = 1 / beta[1, :].transpose() + beta, _, _, _ = np.linalg.lstsq(X, B) + t2s = 1. / beta[1, :].transpose() s0 = np.exp(beta[0, :]).transpose() t2s[np.isinf(t2s)] = 500. s0[np.isnan(s0)] = 0. - t2ss[:, :, :, ne - 2] = np.squeeze(unmask(t2s, mask)) - s0vs[:, :, :, ne - 2] = np.squeeze(unmask(s0, mask)) + t2ss[:, :, :, i_echo-2] = np.squeeze(unmask(t2s, mask)) + s0vs[:, :, :, i_echo-2] = np.squeeze(unmask(s0, mask)) # Limited T2* and S0 maps - fl = np.zeros([nx, ny, nz, len(tes) - 2 + 1]) - for ne in range(Ne - 1): - fl_ = np.squeeze(fl[:, :, :, ne]) - fl_[masksum == ne + 2] = True - fl[:, :, :, ne] = fl_ - fl = np.array(fl, dtype=bool) + fl = np.zeros([nx, ny, nz, len(tes)-1], bool) + for i_echo in range(n_echoes - 1): + fl_ = np.squeeze(fl[:, :, :, i_echo]) + fl_[masksum == i_echo + 2] = True + fl[:, :, :, i_echo] = fl_ t2sa = np.squeeze(unmask(t2ss[fl], masksum > 1)) s0va = np.squeeze(unmask(s0vs[fl], masksum > 1)) @@ -63,39 +133,60 @@ def t2sadmap(catd, mask, tes, masksum, start_echo): def optcom(data, t2, tes, mask, combmode, useG=False): """ - out = optcom(data,t2s) - + Optimally combine BOLD data across TEs. - Input: - - data.shape = (nx,ny,nz,Ne,Nt) - t2s.shape = (nx,ny,nz) - tes.shape = len(Ne) - - Output: + out = optcom(data,t2s) - out.shape = (nx,ny,nz,Nt) + Parameters + ---------- + data : :obj:`numpy.ndarray` + Concatenated BOLD data. Has shape (nx, ny, nz, n_echoes, n_trs) + t2 : :obj:`numpy.ndarray` + 3D map of estimated T2* values. Has shape (nx, ny, nz) + tes : :obj:`numpy.ndarray` + Array of TEs, in seconds. + mask : :obj:`numpy.ndarray` + Brain mask in 3D array. Has shape (nx, ny, nz) + combmode : :obj:`str` + How to combine data. Either 'ste' or 't2s'. + useG : :obj:`bool`, optional + Use G. Default is False. + + Returns + ------- + out : :obj:`numpy.ndarray` + Optimally combined data. Has shape (nx, ny, nz, n_trs) """ - nx, ny, nz, Ne, Nt = data.shape + _, _, _, _, n_trs = data.shape if useG: fdat = fmask(data, mask) ft2s = fmask(t2, mask) - else: fdat = fmask(data, mask) ft2s = fmask(t2, mask) tes = np.array(tes) tes = tes[np.newaxis, :] - ft2s = ft2s[:, np.newaxis] + + if len(t2.shape) == 3: + print('Optimally combining with voxel-wise T2 estimates') + ft2s = ft2s[:, np.newaxis] + else: + print('Optimally combining with voxel- and volume-wise T2 estimates') + ft2s = ft2s[:, :, np.newaxis] if combmode == 'ste': alpha = fdat.mean(-1) * tes else: alpha = tes * np.exp(-tes / ft2s) - alpha = np.tile(alpha[:, :, np.newaxis], (1, 1, Nt)) + if len(t2.shape) == 3: + alpha = np.tile(alpha[:, :, np.newaxis], (1, 1, n_trs)) + else: + alpha = np.swapaxes(alpha, 1, 2) + ax0_idx, ax2_idx = np.where(np.all(alpha == 0, axis=1)) + alpha[ax0_idx, :, ax2_idx] = 1. fout = np.average(fdat, axis=1, weights=alpha) out = unmask(fout, mask) @@ -104,6 +195,14 @@ def optcom(data, t2, tes, mask, combmode, useG=False): def main(options): """ + Estimate T2 and S0, and optimally combine data across TEs. + + Parameters + ---------- + options + label + tes + data """ if options.label is not None: suf = '_%s' % str(options.label) @@ -111,30 +210,26 @@ def main(options): suf = '' tes = [float(te) for te in options.tes] - ne = len(tes) + n_echoes = len(tes) catim = nib.load(options.data[0]) head = catim.get_header() head.extensions = [] head.set_sform(head.get_sform(), code=1) aff = catim.get_affine() - catd = cat2echos(catim.get_data(), ne) - nx, ny, nz, Ne, nt = catd.shape + catd = cat2echos(catim.get_data(), n_echoes) + nx, ny, nz, n_echoes, n_trs = catd.shape print("++ Computing Mask") - mask, masksum = makeadmask(catd, min=False, getsum=True) + mask, masksum = makeadmask(catd, minimum=False, getsum=True) + niwrite(masksum, aff, 'masksum%s.nii' % suf) print("++ Computing Adaptive T2* map") t2s, s0, t2ss, s0vs, t2saf, s0vaf = t2sadmap(catd, mask, tes, masksum, 2) - niwrite(masksum, aff, 'masksum%s.nii' % suf) niwrite(t2ss, aff, 't2ss%s.nii' % suf) niwrite(s0vs, aff, 's0vs%s.nii' % suf) print("++ Computing optimal combination") - tsoc = np.array(optcom(catd, - t2s, - tes, - mask, - options.combmode), + tsoc = np.array(optcom(catd, t2s, tes, mask, options.combmode), dtype=float) # Clean up numerical errors