Skip to content

Commit

Permalink
[REF] Update t2smap workflow for fMRIPrep compatibility (#557)
Browse files Browse the repository at this point in the history 
* Add out-dir argument to t2smap workflow.

Also move logging setup from _main into t2smap_workflow.

* Change t2smap workflow output files.

* I hate tests.

* Oh my tests are the worst.

* Fix the CLI test as well.
  • Loading branch information
@tsalo
tsalo authored Apr 25, 2020
1 parent a2de154 commit cb11c9d
Show file tree
Hide file tree
Showing 2 changed files with 103 additions and 114 deletions.
78 changes: 39 additions & 39 deletions tedana/tests/test_t2smap.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,21 +21,21 @@ def test_basic_t2smap1(self):
data = [op.join(data_dir, 'echo1.nii.gz'),
op.join(data_dir, 'echo2.nii.gz'),
op.join(data_dir, 'echo3.nii.gz')]
workflows.t2smap_workflow(data, [14.5, 38.5, 62.5], combmode='t2s',
fitmode='all', label='t2smap')
out_dir = 'TED.echo1.t2smap'
workflows.t2smap_workflow(data, [14.5, 38.5, 62.5], combmode='t2s',
fitmode='all', out_dir=out_dir)

# Check outputs
assert op.isfile(op.join(out_dir, 'ts_OC.nii.gz'))
img = nib.load(op.join(out_dir, 't2sv.nii.gz'))
assert op.isfile(op.join(out_dir, 'desc-optcom_bold.nii.gz'))
img = nib.load(op.join(out_dir, 'T2StarMap.nii.gz'))
assert len(img.shape) == 3
img = nib.load(op.join(out_dir, 's0v.nii.gz'))
img = nib.load(op.join(out_dir, 'S0Map.nii.gz'))
assert len(img.shape) == 3
img = nib.load(op.join(out_dir, 't2svG.nii.gz'))
img = nib.load(op.join(out_dir, 'desc-full_T2StarMap.nii.gz'))
assert len(img.shape) == 3
img = nib.load(op.join(out_dir, 's0vG.nii.gz'))
img = nib.load(op.join(out_dir, 'desc-full_S0Map.nii.gz'))
assert len(img.shape) == 3
img = nib.load(op.join(out_dir, 'ts_OC.nii.gz'))
img = nib.load(op.join(out_dir, 'desc-optcom_bold.nii.gz'))
assert len(img.shape) == 4

def test_basic_t2smap2(self):
Expand All @@ -47,21 +47,21 @@ def test_basic_t2smap2(self):
data = [op.join(data_dir, 'echo1.nii.gz'),
op.join(data_dir, 'echo2.nii.gz'),
op.join(data_dir, 'echo3.nii.gz')]
workflows.t2smap_workflow(data, [14.5, 38.5, 62.5], combmode='t2s',
fitmode='ts', label='t2smap')
out_dir = 'TED.echo1.t2smap'
workflows.t2smap_workflow(data, [14.5, 38.5, 62.5], combmode='t2s',
fitmode='ts', out_dir=out_dir)

# Check outputs
assert op.isfile(op.join(out_dir, 'ts_OC.nii.gz'))
img = nib.load(op.join(out_dir, 't2sv.nii.gz'))
assert op.isfile(op.join(out_dir, 'desc-optcom_bold.nii.gz'))
img = nib.load(op.join(out_dir, 'T2StarMap.nii.gz'))
assert len(img.shape) == 4
img = nib.load(op.join(out_dir, 's0v.nii.gz'))
img = nib.load(op.join(out_dir, 'S0Map.nii.gz'))
assert len(img.shape) == 4
img = nib.load(op.join(out_dir, 't2svG.nii.gz'))
img = nib.load(op.join(out_dir, 'desc-full_T2StarMap.nii.gz'))
assert len(img.shape) == 4
img = nib.load(op.join(out_dir, 's0vG.nii.gz'))
img = nib.load(op.join(out_dir, 'desc-full_S0Map.nii.gz'))
assert len(img.shape) == 4
img = nib.load(op.join(out_dir, 'ts_OC.nii.gz'))
img = nib.load(op.join(out_dir, 'desc-optcom_bold.nii.gz'))
assert len(img.shape) == 4

def test_basic_t2smap3(self):
Expand All @@ -73,21 +73,21 @@ def test_basic_t2smap3(self):
data = [op.join(data_dir, 'echo1.nii.gz'),
op.join(data_dir, 'echo2.nii.gz'),
op.join(data_dir, 'echo3.nii.gz')]
workflows.t2smap_workflow(data, [14.5, 38.5, 62.5], combmode='paid',
fitmode='all', label='t2smap')
out_dir = 'TED.echo1.t2smap'
workflows.t2smap_workflow(data, [14.5, 38.5, 62.5], combmode='paid',
fitmode='all', out_dir=out_dir)

# Check outputs
assert op.isfile(op.join(out_dir, 'ts_OC.nii.gz'))
img = nib.load(op.join(out_dir, 't2sv.nii.gz'))
assert op.isfile(op.join(out_dir, 'desc-optcom_bold.nii.gz'))
img = nib.load(op.join(out_dir, 'T2StarMap.nii.gz'))
assert len(img.shape) == 3
img = nib.load(op.join(out_dir, 's0v.nii.gz'))
img = nib.load(op.join(out_dir, 'S0Map.nii.gz'))
assert len(img.shape) == 3
img = nib.load(op.join(out_dir, 't2svG.nii.gz'))
img = nib.load(op.join(out_dir, 'desc-full_T2StarMap.nii.gz'))
assert len(img.shape) == 3
img = nib.load(op.join(out_dir, 's0vG.nii.gz'))
img = nib.load(op.join(out_dir, 'desc-full_S0Map.nii.gz'))
assert len(img.shape) == 3
img = nib.load(op.join(out_dir, 'ts_OC.nii.gz'))
img = nib.load(op.join(out_dir, 'desc-optcom_bold.nii.gz'))
assert len(img.shape) == 4

def test_basic_t2smap4(self):
Expand All @@ -99,21 +99,21 @@ def test_basic_t2smap4(self):
data = [op.join(data_dir, 'echo1.nii.gz'),
op.join(data_dir, 'echo2.nii.gz'),
op.join(data_dir, 'echo3.nii.gz')]
workflows.t2smap_workflow(data, [14.5, 38.5, 62.5], combmode='paid',
fitmode='ts', label='t2smap')
out_dir = 'TED.echo1.t2smap'
workflows.t2smap_workflow(data, [14.5, 38.5, 62.5], combmode='paid',
fitmode='ts', out_dir=out_dir)

# Check outputs
assert op.isfile(op.join(out_dir, 'ts_OC.nii.gz'))
img = nib.load(op.join(out_dir, 't2sv.nii.gz'))
assert op.isfile(op.join(out_dir, 'desc-optcom_bold.nii.gz'))
img = nib.load(op.join(out_dir, 'T2StarMap.nii.gz'))
assert len(img.shape) == 4
img = nib.load(op.join(out_dir, 's0v.nii.gz'))
img = nib.load(op.join(out_dir, 'S0Map.nii.gz'))
assert len(img.shape) == 4
img = nib.load(op.join(out_dir, 't2svG.nii.gz'))
img = nib.load(op.join(out_dir, 'desc-full_T2StarMap.nii.gz'))
assert len(img.shape) == 4
img = nib.load(op.join(out_dir, 's0vG.nii.gz'))
img = nib.load(op.join(out_dir, 'desc-full_S0Map.nii.gz'))
assert len(img.shape) == 4
img = nib.load(op.join(out_dir, 'ts_OC.nii.gz'))
img = nib.load(op.join(out_dir, 'desc-optcom_bold.nii.gz'))
assert len(img.shape) == 4

def test_t2smap_cli(self):
Expand All @@ -124,22 +124,22 @@ def test_t2smap_cli(self):
data = [op.join(data_dir, 'echo1.nii.gz'),
op.join(data_dir, 'echo2.nii.gz'),
op.join(data_dir, 'echo3.nii.gz')]
out_dir = 'TED.echo1.t2smap'
args = (['-d'] + data +
['-e', '14.5', '38.5', '62.5', '--combmode', 't2s',
'--fitmode', 'all', '--label', 't2smap'])
'--fitmode', 'all', '--out-dir', out_dir])
workflows.t2smap._main(args)
out_dir = 'TED.echo1.t2smap'

# Check outputs
img = nib.load(op.join(out_dir, 't2sv.nii.gz'))
img = nib.load(op.join(out_dir, 'T2StarMap.nii.gz'))
assert len(img.shape) == 3
img = nib.load(op.join(out_dir, 's0v.nii.gz'))
img = nib.load(op.join(out_dir, 'S0Map.nii.gz'))
assert len(img.shape) == 3
img = nib.load(op.join(out_dir, 't2svG.nii.gz'))
img = nib.load(op.join(out_dir, 'desc-full_T2StarMap.nii.gz'))
assert len(img.shape) == 3
img = nib.load(op.join(out_dir, 's0vG.nii.gz'))
img = nib.load(op.join(out_dir, 'desc-full_S0Map.nii.gz'))
assert len(img.shape) == 3
img = nib.load(op.join(out_dir, 'ts_OC.nii.gz'))
img = nib.load(op.join(out_dir, 'desc-optcom_bold.nii.gz'))
assert len(img.shape) == 4

def teardown_method(self):
Expand Down
139 changes: 64 additions & 75 deletions tedana/workflows/t2smap.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,12 @@ def _get_parser():
type=float,
help='Echo times (in ms). E.g., 15.0 39.0 63.0',
required=True)
optional.add_argument('--out-dir',
dest='out_dir',
type=str,
metavar='PATH',
help='Output directory.',
default='.')
optional.add_argument('--mask',
dest='mask',
metavar='FILE',
Expand All @@ -57,6 +63,17 @@ def _get_parser():
'Dependent ANAlysis. Must be in the same '
'space as `data`.'),
default=None)
optional.add_argument('--fittype',
dest='fittype',
action='store',
choices=['loglin', 'curvefit'],
help='Desired Fitting Method'
'"loglin" means that a linear model is fit'
' to the log of the data, default'
'"curvefit" means that a more computationally'
'demanding monoexponential model is fit'
'to the raw data',
default='loglin')
optional.add_argument('--fitmode',
dest='fitmode',
action='store',
Expand All @@ -74,22 +91,6 @@ def _get_parser():
help=('Combination scheme for TEs: '
't2s (Posse 1999, default), paid (Poser)'),
default='t2s')
optional.add_argument('--label',
dest='label',
type=str,
help='Label for output directory.',
default=None)
optional.add_argument('--fittype',
dest='fittype',
action='store',
choices=['loglin', 'curvefit'],
help='Desired Fitting Method'
'"loglin" means that a linear model is fit'
' to the log of the data, default'
'"curvefit" means that a more computationally'
'demanding monoexponential model is fit'
'to the raw data',
default='loglin')
optional.add_argument('--n-threads',
dest='n_threads',
type=int,
Expand All @@ -114,8 +115,9 @@ def _get_parser():
return parser


def t2smap_workflow(data, tes, mask=None, fitmode='all', combmode='t2s',
label=None, debug=False, fittype='loglin', quiet=False):
def t2smap_workflow(data, tes, out_dir='.', mask=None,
fittype='loglin', fitmode='all', combmode='t2s',
debug=False, quiet=False):
"""
Estimate T2 and S0, and optimally combine data across TEs.
Expand All @@ -126,24 +128,24 @@ def t2smap_workflow(data, tes, mask=None, fitmode='all', combmode='t2s',
list of echo-specific files, in ascending order.
tes : :obj:`list`
List of echo times associated with data in milliseconds.
out_dir : :obj:`str`, optional
Output directory.
mask : :obj:`str`, optional
Binary mask of voxels to include in TE Dependent ANAlysis. Must be spatially
aligned with `data`.
fittype : {'loglin', 'curvefit'}, optional
Monoexponential fitting method.
'loglin' means to use the the default linear fit to the log of
the data.
'curvefit' means to use a monoexponential fit to the raw data,
which is slightly slower but may be more accurate.
fitmode : {'all', 'ts'}, optional
Monoexponential model fitting scheme.
'all' means that the model is fit, per voxel, across all timepoints.
'ts' means that the model is fit, per voxel and per timepoint.
Default is 'all'.
combmode : {'t2s', 'paid'}, optional
Combination scheme for TEs: 't2s' (Posse 1999, default), 'paid' (Poser).
label : :obj:`str` or :obj:`None`, optional
Label for output directory. Default is None.
fittype : {'loglin', 'curvefit'}, optional
Monoexponential fitting method.
'loglin' means to use the the default linear fit to the log of
the data.
'curvefit' means to use a monoexponential fit to the raw data,
which is slightly slower but may be more accurate.
Other Parameters
----------------
Expand All @@ -155,29 +157,37 @@ def t2smap_workflow(data, tes, mask=None, fitmode='all', combmode='t2s',
Notes
-----
This workflow writes out several files, which are written out to a folder
named TED.[ref_label].[label] if ``label`` is provided and TED.[ref_label]
if not. ``ref_label`` is determined based on the name of the first ``data``
file.
Files are listed below:
This workflow writes out several files, which are described below:
====================== =================================================
Filename Content
====================== =================================================
t2sv.nii Limited estimated T2* 3D map or 4D timeseries.
Will be a 3D map if ``fitmode`` is 'all' and a
4D timeseries if it is 'ts'.
s0v.nii Limited S0 3D map or 4D timeseries.
t2svG.nii Full T2* map/timeseries. The difference between
the limited and full maps is that, for voxels
affected by dropout where only one echo contains
good data, the full map uses the single echo's
value while the limited map has a NaN.
s0vG.nii Full S0 map/timeseries.
ts_OC.nii Optimally combined timeseries.
====================== =================================================
========================== =================================================
Filename Content
========================== =================================================
T2StarMap.nii.gz Limited estimated T2* 3D map or 4D timeseries.
Will be a 3D map if ``fitmode`` is 'all' and a
4D timeseries if it is 'ts'.
S0Map.nii.gz Limited S0 3D map or 4D timeseries.
desc-full_T2StarMap.nii.gz Full T2* map/timeseries. The difference between
the limited and full maps is that, for voxels
affected by dropout where only one echo contains
good data, the full map uses the single echo's
value while the limited map has a NaN.
desc-full_S0Map.nii.gz Full S0 map/timeseries.
desc-optcom_bold.nii.gz Optimally combined timeseries.
========================== =================================================
"""
out_dir = op.abspath(out_dir)
if not op.isdir(out_dir):
os.mkdir(out_dir)

if debug and not quiet:
logging.basicConfig(level=logging.DEBUG)
elif quiet:
logging.basicConfig(level=logging.WARNING)
else:
logging.basicConfig(level=logging.INFO)

LGR.info('Using output directory: {}'.format(out_dir))

# ensure tes are in appropriate format
tes = [float(te) for te in tes]
n_echos = len(tes)
Expand All @@ -191,22 +201,6 @@ def t2smap_workflow(data, tes, mask=None, fitmode='all', combmode='t2s',
n_samp, n_echos, n_vols = catd.shape
LGR.debug('Resulting data shape: {}'.format(catd.shape))

try:
ref_label = op.basename(ref_img).split('.')[0]
except (TypeError, AttributeError):
ref_label = op.basename(str(data[0])).split('.')[0]

if label is not None:
out_dir = 'TED.{0}.{1}'.format(ref_label, label)
else:
out_dir = 'TED.{0}'.format(ref_label)
out_dir = op.abspath(out_dir)
if not op.isdir(out_dir):
LGR.info('Creating output directory: {}'.format(out_dir))
os.mkdir(out_dir)
else:
LGR.info('Using output directory: {}'.format(out_dir))

if mask is None:
LGR.info('Computing adaptive mask')
else:
Expand Down Expand Up @@ -243,23 +237,18 @@ def t2smap_workflow(data, tes, mask=None, fitmode='all', combmode='t2s',
s0_limited[s0_limited < 0] = 0
t2s_limited[t2s_limited < 0] = 0

io.filewrite(utils.millisec2sec(t2s_limited), op.join(out_dir, 't2sv.nii'), ref_img)
io.filewrite(s0_limited, op.join(out_dir, 's0v.nii'), ref_img)
io.filewrite(utils.millisec2sec(t2s_full), op.join(out_dir, 't2svG.nii'), ref_img)
io.filewrite(s0_full, op.join(out_dir, 's0vG.nii'), ref_img)
io.filewrite(OCcatd, op.join(out_dir, 'ts_OC.nii'), ref_img)
io.filewrite(utils.millisec2sec(t2s_limited),
op.join(out_dir, 'T2StarMap.nii.gz'), ref_img)
io.filewrite(s0_limited, op.join(out_dir, 'S0Map.nii.gz'), ref_img)
io.filewrite(utils.millisec2sec(t2s_full),
op.join(out_dir, 'desc-full_T2StarMap.nii.gz'), ref_img)
io.filewrite(s0_full, op.join(out_dir, 'desc-full_S0Map.nii.gz'), ref_img)
io.filewrite(OCcatd, op.join(out_dir, 'desc-optcom_bold.nii.gz'), ref_img)


def _main(argv=None):
"""T2smap entry point"""
options = _get_parser().parse_args(argv)
if options.debug and not options.quiet:
logging.basicConfig(level=logging.DEBUG)
elif options.quiet:
logging.basicConfig(level=logging.WARNING)
else:
logging.basicConfig(level=logging.INFO)

kwargs = vars(options)
n_threads = kwargs.pop('n_threads')
n_threads = None if n_threads == -1 else n_threads
Expand Down

0 comments on commit cb11c9d

Please sign in to comment.