Skip to content

Commit

Permalink
Merge pull request #1 from TomMaullin/refactorFit
Browse files Browse the repository at this point in the history 
Refactor Fit
  • Loading branch information
@tsalo
tsalo authored May 14, 2018
2 parents 8680466 + d887611 commit 6cc8f49
Showing 1 changed file with 59 additions and 75 deletions.
134 changes: 59 additions & 75 deletions tedana/model/monoexponential.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@
LGR = logging.getLogger(__name__)


def fit_decay(data, tes, mask, masksum, start_echo):
def fit_decay(data, tes, mask, masksum, start_echo, dim):
"""
Fit voxel-wise monoexponential decay models to estimate T2* and S0 maps.
Expand All @@ -30,6 +30,8 @@ def fit_decay(data, tes, mask, masksum, start_echo):
given sample
start_echo : int
First echo to consider
dim : The dimensions we wish to consider - e.g. if dim is 4 we only loop over
x by y by z by E. If dim is 5 we loop over x by y by z by E by T.
Returns
-------
Expand Down Expand Up @@ -58,39 +60,56 @@ def fit_decay(data, tes, mask, masksum, start_echo):
3. Generate limited :math:`T_2^*` and :math:`S_0` maps by doing something.
"""

n_samp, n_echos, n_vols = data.shape
data = data[mask]
t2ss, s0vs = np.zeros([n_samp, n_echos - 1]), np.zeros([n_samp, n_echos - 1])

for echo in range(start_echo, n_echos + 1):
# perform log linear fit of echo times against MR signal
# make DV matrix: samples x (time series * echos)
B = np.log((np.abs(data[:, :echo, :]) + 1).reshape(len(data), -1).T)
# make IV matrix: intercept/TEs x (time series * echos)
x = np.column_stack([np.ones(echo), [-te for te in tes[:echo]]])
X = np.repeat(x, n_vols, axis=0)

beta = np.linalg.lstsq(X, B)[0]
t2s = 1. / beta[1, :].T
s0 = np.exp(beta[0, :]).T

t2s[np.isinf(t2s)] = 500. # why 500?
s0[np.isnan(s0)] = 0. # why 0?

t2ss[..., echo - 2] = np.squeeze(utils.unmask(t2s, mask))
s0vs[..., echo - 2] = np.squeeze(utils.unmask(s0, mask))

# create limited T2* and S0 maps
fl = np.zeros([n_samp, len(tes) - 1], dtype=bool)
for echo in range(n_echos - 1):
fl_ = np.squeeze(fl[..., echo])
fl_[masksum == echo + 2] = True
fl[..., echo] = fl_
t2sa, s0va = utils.unmask(t2ss[fl], masksum > 1), utils.unmask(s0vs[fl], masksum > 1)
# t2sa[masksum > 1], s0va[masksum > 1] = t2ss[fl], s0vs[fl]

# create full T2* maps with S0 estimation errors
t2saf, s0vaf = t2sa.copy(), s0va.copy()
if dim == 5:
nx, ny, nz, n_echoes, n_trs = data.shape
else:
nx, ny, nz, n_echoes,_ = data.shape
n_trs = 1

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()

# 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[:, :i_echo, :]) + 1,
(n_voxels, i_echo*n_trs)).transpose()
B = np.log(B)
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, _, _, _ = 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]

Expand Down Expand Up @@ -133,8 +152,8 @@ def fit_decay_ts(data, mask, tes, masksum, start_echo):
s0vaf : (S x E x T) :obj:`numpy.ndarray`
Full S0 map
"""
nx, ny, nz, n_echos, n_trs = data.shape
echodata = data[mask]
nx, ny, nz, n_echoes, n_trs = data.shape
echodata = fmask(data, mask)
tes = np.array(tes)

t2sa_ts = np.zeros([nx, ny, nz, n_trs])
Expand All @@ -143,48 +162,13 @@ def fit_decay_ts(data, mask, tes, masksum, start_echo):
s0vaf_ts = np.zeros([nx, ny, nz, n_trs])

for vol in range(echodata.shape[-1]):
t2ss = np.zeros([nx, ny, nz, n_echos - 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_echos + 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(utils.unmask(t2s, mask))
s0vs[:, :, :, i_echo-2] = np.squeeze(utils.unmask(s0, mask))

# Limited T2* and S0 maps
fl = np.zeros([nx, ny, nz, len(tes)-1], bool)
for i_echo in range(n_echos - 1):
fl_ = np.squeeze(fl[:, :, :, i_echo])
fl_[masksum == i_echo + 2] = True
fl[:, :, :, i_echo] = fl_
t2sa = np.squeeze(utils.unmask(t2ss[fl], masksum > 1))
s0va = np.squeeze(utils.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, s0va, t2ss, s0vs, t2saf, s0vaf] = fit_decay(
data, mask, tes, masksum, start_echo, 4)

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
return t2sa_ts, s0va_ts, t2saf_ts, s0vaf_ts

0 comments on commit 6cc8f49

Please sign in to comment.