[ENH] Adds extraction of physio signals from DICOMs

Dockerfile

@@ -62,6 +62,12 @@ RUN apt-get update -qq \
     && make install \
     && rm -rf /tmp/dcm2niix
 
+# Download bidsphysio (we'll install later):
+RUN mkdir -p /src/bidsphysio \
+    && curl -sSL https://github.com/cbinyu/bidsphysio/archive/v3.0.tar.gz \
+        | tar -vxz -C /src/bidsphysio --exclude='bidsphysio/tests'\    
+ --strip-components=1
+
 RUN apt-get update -qq \
     && apt-get install -y -q --no-install-recommends \
            git \
@@ -99,7 +105,8 @@ RUN export PATH="/opt/miniconda-latest/bin:$PATH" \
     && sync && conda clean -tipsy && sync \
     && bash -c "source activate base \
     &&   pip install --no-cache-dir --editable \
-             '/src/heudiconv[all]'" \
+             '/src/heudiconv[all]' \
+	     '/src/bidsphysio'" \
     && rm -rf ~/.cache/pip/* \
     && sync
 

heudiconv/convert.py

@@ -288,6 +288,8 @@ def convert(items, converter, scaninfo_suffix, custom_callable, with_prov,
             if outtype == 'dicom':
                 convert_dicom(item_dicoms, bids_options, prefix,
                               outdir, tempdirs, symlink, overwrite)
+            elif outtype == 'physio':
+                convert_physio(item_dicoms, bids_options, prefix)
             elif outtype in ['nii', 'nii.gz']:
                 assert converter == 'dcm2niix', ('Invalid converter '
                                                  '{}'.format(converter))
@@ -407,6 +409,43 @@ def convert_dicom(item_dicoms, bids_options, prefix,
                 shutil.copyfile(filename, outfile)
 
 
+def convert_physio(item_dicoms, bids_options, prefix):
+    """Save DICOM physiology as BIDS physio files
+    Parameters
+    ----------
+    item_dicoms : list of filenames
+        DICOMs to save
+    bids_options : list or None
+        If not None then save to BIDS format. List may be empty
+        or contain bids specific options
+    prefix : string
+        Conversion outname
+    Returns
+    -------
+    None
+    """
+    if bids_options is not None:
+        try:
+            from bidsphysio import dcm2bidsphysio
+        except ImportError:
+            lgr.warning(
+                "Dcm2bidsphysio not found. "
+                "Not extracting physiological recordings."
+            )
+            return
+
+        item_dicoms = list(map(op.abspath, item_dicoms)) # absolute paths
+        if len(item_dicoms)>1:
+            lgr.warning(
+                "More than one PHYSIO file has been found for this series. "
+                "If each file corresponds to a different signal, all is OK. "
+                "If multiple files have the same signal, only the signal "
+                "from the last file will be saved."
+            )
+        for dicom_file in item_dicoms:
+            dcm2bidsphysio.dcm2bids( dicom_file, prefix )
+
+
 def nipype_convert(item_dicoms, prefix, with_prov, bids_options, tmpdir, dcmconfig=None):
     """
     Converts DICOMs grouped from heuristic using Nipype's Dcm2niix interface.

heudiconv/dicoms.py

@@ -265,8 +265,13 @@ def group_dicoms_into_seqinfos(files, grouping, file_filter=None,
         series_id = '-'.join(map(str, series_id))
         if mw.image_shape is None:
             # this whole thing has no image data (maybe just PSg DICOMs)
-            # nothing to see here, just move on
-            continue
+            # If this is a Siemens PhysioLog, keep it:
+            if mw.dcm_data.SeriesDescription.endswith('_PhysioLog'):
+                # just give it a dummy shape, so that we can continue:
+                mw.image_shape=(0,0,0)
+            else:
+                # nothing to see here, just move on
+                continue
         seqinfo = create_seqinfo(mw, series_files, series_id)
 
         if per_studyUID:

heudiconv/heuristics/bids_physio.py

@@ -0,0 +1,99 @@
+"""
+Heuristic demonstrating extraction of physiological data from CMRR
+fMRI DICOMs
+
+We want to make sure the run number for the _sbref, _phase and 
+_physio matches that of the corresponding _bold. For "normal"
+scanning, you can just rely on the {item} value, but if you have a
+functional run with just saving the magnitude and then one saving
+both magnitude and phase, you would have _run-01_bold, _run-02_bold
+and _run-01_phase, but the phase image corresponds to _run-02_bold,
+so the run number in the filename will not match
+"""
+
+
+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
+    """
+
+    info = {}
+    run_no = 0
+    for idx, s in enumerate(seqinfo):
+        # We want to make sure the _SBRef, PhysioLog and phase series
+        # (if present) are labeled the same as the main (magnitude)
+        # image. So we only focus on the magnitude series (to exclude
+        # phase images) without _SBRef at the end of the series_
+        # description and then we search if the phase and/or _SBRef
+        # are present.
+        if (
+            'epfid2d' in s.sequence_name
+            and (
+                'M' in s.image_type
+                or 'FMRI' in s.image_type
+            )
+            and not s.series_description.lower().endswith('_sbref')
+            and not 'DERIVED' in s.image_type
+        ):
+            run_no += 1
+            bold = create_key(
+                'sub-{subject}/func/sub-{subject}_task-test_run-%02d_bold' % run_no
+            )
+            info[bold] = [{'item': s.series_id}]
+            next_series = idx+1    # used for physio log below
+
+            ###   is phase image present?   ###
+            # At least for Siemens systems, if magnitude/phase was
+            # selected, the phase images come as a separate series
+            # immediatelly following the magnitude series.
+            # (note: make sure you don't check beyond the number of
+            # elements in seqinfo...)
+            if (
+                idx+1 < len(seqinfo)
+                and 'P' in seqinfo[idx+1].image_type
+            ):
+                phase = create_key(
+                    'sub-{subject}/func/sub-{subject}_task-test_run-%02d_phase' % run_no
+                )
+                info[phase] = [{'item': seqinfo[idx+1].series_id}]
+                next_series = idx+2    # used for physio log below
+
+            ###   SBREF   ###
+            # here, within the functional run code, check to see if
+            # the previous run's series_description ended in _sbref,
+            # to assign the same run number.
+            if (
+                idx > 0
+                and seqinfo[idx-1].series_description.lower().endswith('_sbref')
+            ):
+                sbref = create_key(
+                    'sub-{subject}/func/sub-{subject}_task-test_run-%02d_sbref' % run_no
+                )
+                info[sbref] = [{'item': seqinfo[idx-1].series_id}]
+
+            ###   PHYSIO LOG   ###
+            # here, within the functional run code, check to see if
+            # the next run image_type lists "PHYSIO", to assign the
+            # same run number.
+            if (
+                next_series < len(seqinfo)
+                and 'PHYSIO' in seqinfo[next_series].image_type
+            ):
+                physio = create_key(
+                    'sub-{subject}/func/sub-{subject}_task-test_run-%02d_physio' % run_no,
+                    outtype = ('physio',)
+                )
+                info[physio] = [{'item': seqinfo[next_series].series_id}]
+
+    return info

heudiconv/tests/data/samplePhysio/README.txt

@@ -0,0 +1,4 @@
+samplePhysio dataset
+
+It contains phantom functional images and physiological recordings
+using CMRR Multi-Band EPI saved as DICOMs.

heudiconv/tests/test_regression.py

@@ -17,6 +17,12 @@
 except ImportError:
     have_datalad = False
 
+have_bidsphysio = True
+try:
+    from bidsphysio import dcm2bidsphysio
+except ImportError:
+    have_bidsphysio = False
+
 
 @pytest.mark.skipif(not have_datalad, reason="no datalad")
 @pytest.mark.parametrize('subject', ['sub-sid000143'])
@@ -84,6 +90,32 @@ def test_multiecho(tmpdir, subject='MEEPI', heuristic='bids_ME.py'):
         assert 'echo-' not in event
 
 
+@pytest.mark.skipif(not have_bidsphysio, reason="no bidsphysio")
+def test_physio(tmpdir, subject='samplePhysio', heuristic='bids_physio.py'):
+    tmpdir.chdir()
+    outdir = tmpdir.mkdir('out').strpath
+    template = op.join("{subject}/*/*.dcm")
+    args = gen_heudiconv_args(
+        TESTS_DATA_PATH, outdir, subject, heuristic,template=template
+    )
+    runner(args)  # run conversion
+
+    # Check we get only one image file:
+    func_images = glob(op.join('out', 'sub-' + subject, 'func', '*.nii.gz'))
+    assert len(func_images) == 1
+    # The corresponding json:
+    _json = func_images[0].replace('.nii.gz', '.json')
+    assert op.exists(_json)
+    # For each physiological signal, we get the json and tsv.gz:
+    for s in ['respiratory','cardiac']:
+        expectedFileName = func_images[0].replace(
+            '_bold.nii.gz',
+            '_recording-' + s + '_physio'
+        )
+        assert op.exists(expectedFileName + '.json')
+        assert op.exists(expectedFileName + '.tsv.gz')
+
+
 @pytest.mark.parametrize('subject', ['merged'])
 def test_grouping(tmpdir, subject):
     dicoms = [