Merge pull request #293 from mgxd/multiecho

[ENH,RF] Multiecho support redux & --files usage changes

Author: yarikoptic

heudiconv/bids.py

@@ -25,6 +25,11 @@
 
 lgr = logging.getLogger(__name__)
 
+
+class BIDSError(Exception):
+    pass
+
+
 def populate_bids_templates(path, defaults={}):
     """Premake BIDS text files with templates"""
 
@@ -113,6 +118,25 @@ def populate_aggregated_jsons(path):
         # create a stub onsets file for each one of those
         suf = '_bold.json'
         assert fpath.endswith(suf)
+        # specify the name of the '_events.tsv' file:
+        if '_echo-' in fpath:
+            # multi-echo sequence: bids (1.1.0) specifies just one '_events.tsv'
+            #   file, common for all echoes.  The name will not include _echo-.
+            # TODO: RF to use re.match for better readability/robustness
+            # So, find out the echo number:
+            fpath_split = fpath.split('_echo-', 1)         # split fpath using '_echo-'
+            fpath_split_2 = fpath_split[1].split('_', 1)   # split the second part of fpath_split using '_'
+            echoNo = fpath_split_2[0]                      # get echo number
+            if echoNo == '1':
+                if len(fpath_split_2) != 2:
+                    raise ValueError("Found no trailer after _echo-")
+                # we modify fpath to exclude '_echo-' + echoNo:
+                fpath = fpath_split[0] + '_' + fpath_split_2[1]
+            else:
+                # for echoNo greater than 1, don't create the events file, so go to
+                #   the next for loop iteration:
+                continue
+
         events_file = fpath[:-len(suf)] + '_events.tsv'
         # do not touch any existing thing, it may be precious
         if not op.lexists(events_file):

heudiconv/convert.py

@@ -24,6 +24,7 @@
     save_scans_key,
     tuneup_bids_json_files,
     add_participant_record,
+    BIDSError
 )
 from .dicoms import (
     group_dicoms_into_seqinfos,
@@ -243,7 +244,7 @@ def convert(items, converter, scaninfo_suffix, custom_callable, with_prov,
         if not isinstance(outtypes, (list, tuple)):
             outtypes = (outtypes,)
 
-        prefix_dirname = op.dirname(prefix + '.ext')
+        prefix_dirname = op.dirname(prefix)
         outname_bids = prefix + '.json'
         bids_outfiles = []
         lgr.info('Converting %s (%d DICOMs) -> %s . '
@@ -462,6 +463,8 @@ def save_converted_files(res, item_dicoms, bids, outtype, prefix, outname_bids,
     """
     from nipype.interfaces.base import isdefined
 
+    prefix_dirname, prefix_basename = op.split(prefix)
+
     bids_outfiles = []
     res_files = res.outputs.converted_files
 
@@ -492,16 +495,112 @@ def save_converted_files(res, item_dicoms, bids, outtype, prefix, outname_bids,
         # Also copy BIDS files although they might need to
         # be merged/postprocessed later
         bids_files = sorted(res.outputs.bids
-                      if len(res.outputs.bids) == len(res_files)
-                      else [None] * len(res_files))
+                     if len(res.outputs.bids) == len(res_files)
+                     else [None] * len(res_files))
+
+        ###   Do we have a multi-echo series?   ###
+        #   Some Siemens sequences (e.g. CMRR's MB-EPI) set the label 'TE1',
+        #   'TE2', etc. in the 'ImageType' field. However, other seqs do not
+        #   (e.g. MGH ME-MPRAGE). They do set a 'EchoNumber', but not for the
+        #   first echo.  To compound the problem, the echoes are NOT in order,
+        #   so the first NIfTI file does not correspond to echo-1, etc. So, we
+        #   need to know, beforehand, whether we are dealing with a multi-echo
+        #   series. To do that, the most straightforward way is to read the
+        #   echo times for all bids_files and see if they are all the same or not.
+
+        # Check for echotime information
+        echo_times = set()
+
+        for bids_file in bids_files:
+            if bids_file:
+                # check for varying EchoTimes
+                echot = load_json(bids_file).get('EchoTime', None)
+                if echot is not None:
+                    echo_times.add(echot)
 
+        # To see if the echo times are the same, convert it to a set and see if
+        # only one remains:
+        is_multiecho = len(echo_times) >= 1 if echo_times else False
+
+        ### Loop through the bids_files, set the output name and save files
         for fl, suffix, bids_file in zip(res_files, suffixes, bids_files):
-            outname = "%s%s.%s" % (prefix, suffix, outtype)
-            safe_copyfile(fl, outname, overwrite)
+
+            # TODO: monitor conversion duration
             if bids_file:
-                outname_bids_file = "%s%s.json" % (prefix, suffix)
+                fileinfo = load_json(bids_file)
+
+            # set the prefix basename for this specific file (we'll modify it,
+            # and we don't want to modify it for all the bids_files):
+            this_prefix_basename = prefix_basename
+
+            # _sbref sequences reconstructing magnitude and phase generate
+            # two NIfTI files IN THE SAME SERIES, so we cannot just add
+            # the suffix, if we want to be bids compliant:
+            if bids_file and this_prefix_basename.endswith('_sbref'):
+                # Check to see if it is magnitude or phase reconstruction:
+                if 'M' in fileinfo.get('ImageType'):
+                    mag_or_phase = 'magnitude'
+                elif 'P' in fileinfo.get('ImageType'):
+                    mag_or_phase = 'phase'
+                else:
+                    mag_or_phase = suffix
+
+                # Insert reconstruction label
+                if not ("_rec-%s" % mag_or_phase) in this_prefix_basename:
+
+                    # If "_rec-" is specified, prepend the 'mag_or_phase' value.
+                    if ('_rec-' in this_prefix_basename):
+                        raise BIDSError(
+                        "Reconstruction label for multi-echo single-band"
+                        " reference images will be automatically set, remove"
+                        " from heuristic"
+                        )
+
+                    # If not, insert "_rec-" + 'mag_or_phase' into the prefix_basename
+                    #   **before** "_run", "_echo" or "_sbref", whichever appears first:
+                    for label in ['_run', '_echo', '_sbref']:
+                        if (label in this_prefix_basename):
+                            this_prefix_basename = this_prefix_basename.replace(
+                                label, "_rec-%s%s" % (mag_or_phase, label)
+                            )
+                            break
+
+            # Now check if this run is multi-echo
+            # (Note: it can be _sbref and multiecho, so don't use "elif"):
+            # For multi-echo sequences, we have to specify the echo number in
+            # the file name:
+            if bids_file and is_multiecho:
+                # Get the EchoNumber from json file info.  If not present, it's echo-1
+                echo_number = fileinfo.get('EchoNumber', 1)
+
+
+                supported_multiecho = ['_bold', '_epi', '_sbref', '_T1w']
+                # Now, decide where to insert it.
+                # Insert it **before** the following string(s), whichever appears first.
+                for imgtype in supported_multiecho:
+                    if (imgtype in this_prefix_basename):
+                        this_prefix_basename = this_prefix_basename.replace(
+                            imgtype, "_echo-%d%s" % (echo_number, imgtype)
+                        )
+                        break
+
+            # Fallback option:
+            # If we have failed to modify this_prefix_basename, because it didn't fall
+            #   into any of the options above, just add the suffix at the end:
+            if this_prefix_basename == prefix_basename:
+                this_prefix_basename += suffix
+
+            # Finally, form the outname by stitching the directory and outtype:
+            outname = op.join(prefix_dirname, this_prefix_basename)
+            outfile = outname + '.' + outtype
+
+            # Write the files needed:
+            safe_copyfile(fl, outfile, overwrite)
+            if bids_file:
+                outname_bids_file = "%s.json" % (outname)
                 safe_copyfile(bids_file, outname_bids_file, overwrite)
                 bids_outfiles.append(outname_bids_file)
+
     # res_files is not a list
     else:
         outname = "{}.{}".format(prefix, outtype)

heudiconv/heuristics/bids_ME.py

@@ -0,0 +1,29 @@
+"""Heuristic demonstrating conversion of the Multi-Echo sequences.
+
+It only cares about converting sequences which have _ME_ in their
+series_description and outputs to BIDS.
+"""
+
+
+def create_key(template, outtype=('nii.gz',), annotation_classes=None):
+    if template is None or not template:
+        raise ValueError('Template must be a valid format string')
+    return template, outtype, annotation_classes
+
+def infotodict(seqinfo):
+    """Heuristic evaluator for determining which runs belong where
+
+    allowed template fields - follow python string module:
+
+    item: index within category
+    subject: participant id
+    seqitem: run number during scanning
+    subindex: sub index within group
+    """
+    bold = create_key('sub-{subject}/func/sub-{subject}_task-test_run-{item}_bold')
+
+    info = {bold: []}
+    for s in seqinfo:
+        if '_ME_' in s.series_description:
+            info[bold].append(s.series_id)
+    return info

heudiconv/parser.py

@@ -167,11 +167,6 @@ def get_study_sessions(dicom_dir_template, files_opt, heuristic, outdir,
             dcmfilter=getattr(heuristic, 'filter_dicom', None),
             grouping=grouping)
 
-        if not getattr(heuristic, 'infotoids', None):
-            raise NotImplementedError(
-                "For now, if no subj template is provided, requiring "
-                "heuristic to have infotoids")
-
         if sids:
             if not (len(sids) == 1 and len(seqinfo_dict) == 1):
                 raise RuntimeError(
@@ -184,6 +179,24 @@ def get_study_sessions(dicom_dir_template, files_opt, heuristic, outdir,
         else:
             sid = None
 
+        if not getattr(heuristic, 'infotoids', None):
+            # allow bypass with subject override
+            if not sid:
+                raise NotImplementedError("Cannot guarantee subject id - add "
+                                          "`infotoids` to heuristic file or "
+                                          "provide `--subjects` option")
+            lgr.warn("Heuristic is missing an `infotoids` method, assigning "
+                     "empty method and using provided subject id %s."
+                     "Provide `session` and `locator` fields for best results."
+                     , sid)
+            def infotoids(seqinfos, outdir):
+                return {
+                    'locator': None,
+                    'session': None,
+                    'subject': None
+                    }
+            heuristic.infotoids = infotoids
+
         for studyUID, seqinfo in seqinfo_dict.items():
             # so we have a single study, we need to figure out its
             # locator, session, subject

heudiconv/tests/test_heuristics.py

@@ -64,10 +64,11 @@ def test_reproin_largely_smoke(tmpdir, heuristic, invocation):
             runner(args + ['--subjects', 'sub1', 'sub2'])
 
         if heuristic != 'reproin':
-            # none other heuristic has mighty infotoids atm
+            # if subject is not overriden, raise error
             with pytest.raises(NotImplementedError):
                 runner(args)
             return
+
     runner(args)
     ds = Dataset(str(tmpdir))
     assert ds.is_installed()

heudiconv/tests/test_regression.py

@@ -1,6 +1,7 @@
 """Testing conversion with conversion saved on datalad"""
 import json
 from glob import glob
+import os.path as op
 
 import pytest
 
@@ -11,8 +12,8 @@
 except ImportError:
     have_datalad = False
 
-import heudiconv
 from heudiconv.cli.run import main as runner
+from heudiconv.utils import load_json
 # testing utilities
 from .utils import fetch_data, gen_heudiconv_args
 
@@ -24,12 +25,19 @@
 def test_conversion(tmpdir, subject, heuristic, anon_cmd):
     tmpdir.chdir()
     try:
-        datadir = fetch_data(tmpdir.strpath, subject)
+        datadir = fetch_data(tmpdir.strpath,
+                             "dbic/QA",  # path from datalad database root
+                             getpath=op.join('sourcedata', subject))
     except IncompleteResultsError as exc:
         pytest.skip("Failed to fetch test data: %s" % str(exc))
     outdir = tmpdir.mkdir('out').strpath
 
-    args = gen_heudiconv_args(datadir, outdir, subject, heuristic, anon_cmd)
+    args = gen_heudiconv_args(datadir,
+                              outdir,
+                              subject,
+                              heuristic,
+                              anon_cmd,
+                              template=op.join('sourcedata/{subject}/*/*/*.tgz'))
     runner(args) # run conversion
 
     # verify functionals were converted
@@ -38,8 +46,38 @@ def test_conversion(tmpdir, subject, heuristic, anon_cmd):
 
     # compare some json metadata
     json_ = '{}/task-rest_acq-24mm64sl1000tr32te600dyn_bold.json'.format
-    orig, conv = (json.load(open(json_(datadir))),
-                  json.load(open(json_(outdir))))
+    orig, conv = (load_json(json_(datadir)),
+                  load_json(json_(outdir)))
     keys = ['EchoTime', 'MagneticFieldStrength', 'Manufacturer', 'SliceTiming']
     for key in keys:
         assert orig[key] == conv[key]
+
+@pytest.mark.skipif(not have_datalad, reason="no datalad")
+def test_multiecho(tmpdir, subject='MEEPI', heuristic='bids_ME.py'):
+    tmpdir.chdir()
+    try:
+        datadir = fetch_data(tmpdir.strpath, "dicoms/velasco/MEEPI")
+    except IncompleteResultsError as exc:
+        pytest.skip("Failed to fetch test data: %s" % str(exc))
+
+    outdir = tmpdir.mkdir('out').strpath
+    args = gen_heudiconv_args(datadir, outdir, subject, heuristic)
+    runner(args) # run conversion
+
+    # check if we have echo functionals
+    echoes = glob(op.join('out', 'sub-' + subject, 'func', '*echo*nii.gz'))
+    assert len(echoes) == 3
+
+    # check EchoTime of each functional
+    # ET1 < ET2 < ET3
+    prev_echo = 0
+    for echo in sorted(echoes):
+        _json = echo.replace('.nii.gz', '.json')
+        assert _json
+        echotime = load_json(_json).get('EchoTime', None)
+        assert echotime > prev_echo
+        prev_echo = echotime
+
+    events = glob(op.join('out', 'sub-' + subject, 'func', '*events.tsv'))
+    for event in events:
+        assert 'echo-' not in event