[ENH] Adds extraction of physio signals from DICOMs

heudiconv/convert.py

@@ -461,6 +461,8 @@
             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))
@@ -580,6 +582,49 @@
                 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 None:
+        return
+
+    try:
+        from bidsphysio.dcm2bids.dcm2bidsphysio import dcm2bids
+    except ImportError:
+        lgr.warning(
+            "bidsphysio.dcm2bids 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:
+        physio_data = dcm2bids(dicom_file)
+        if physio_data.labels():
+            physio_data.save_to_bids_with_trigger(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

@@ -89,7 +89,8 @@
         patient_age=dcminfo.get('PatientAge'),
         patient_sex=dcminfo.get('PatientSex'),
         date=dcminfo.get('AcquisitionDate'),
-        series_uid=dcminfo.get('SeriesInstanceUID')
+        series_uid=dcminfo.get('SeriesInstanceUID'),
+        time=dcminfo.get('AcquisitionTime'),
     )
     return seqinfo
 
@@ -265,8 +266,16 @@
         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 PhoenixZipReport or PhysioLog, keep it:
+            if (
+                mw.dcm_data.SeriesDescription == 'PhoenixZIPReport'
+                or 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_PhoenixReport.py

@@ -0,0 +1,38 @@
+"""Heuristic demonstrating conversion of the PhoenixZIPReport from Siemens.
+
+It only cares about converting a series with have PhoenixZIPReport in their
+series_description and outputs **only to sourcedata**.
+"""
+
+
+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
+    """
+    sbref = create_key('sub-{subject}/func/sub-{subject}_task-QA_sbref', outtype=('nii.gz', 'dicom',))
+    scout = create_key('sub-{subject}/anat/sub-{subject}_T1w', outtype=('nii.gz', 'dicom',))
+    phoenix_doc = create_key('sub-{subject}/misc/sub-{subject}_phoenix', outtype=('dicom',))
+
+    info = {sbref: [], scout: [], phoenix_doc: []}
+    for s in seqinfo:
+        if (
+            'PhoenixZIPReport' in s.series_description
+            and s.image_type[3] == 'CSA REPORT'
+        ):
+            info[phoenix_doc].append({'item': s.series_id})
+        if 'scout' in s.series_description.lower():
+            info[scout].append({'item': s.series_id})
+
+    return info

heudiconv/heuristics/bids_physio.py

@@ -0,0 +1,70 @@
+"""
+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
+
+            ###   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/info.py

@@ -43,7 +43,10 @@
     'extras': [
         'duecredit',  # optional dependency
     ],  # Requires patched version ATM ['dcmstack'],
-    'datalad': ['datalad >=%s' % MIN_DATALAD_VERSION]
+    'datalad': ['datalad >=%s' % MIN_DATALAD_VERSION],
+    'physio': [
+        'bidsphysio.dcm2bids >=1.4.3; python_version>"3.5"',   # if dicoms with physio need to be converted
+    ]
 }
 
 # Flatten the lists

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_dicoms.py

@@ -1,12 +1,18 @@
 import os.path as op
 import json
+from glob import glob
 
 import pytest
 
 from heudiconv.external.pydicom import dcm
 from heudiconv.cli.run import main as runner
 from heudiconv.convert import nipype_convert
-from heudiconv.dicoms import parse_private_csa_header, embed_dicom_and_nifti_metadata
+from heudiconv.dicoms import (
+    OrderedDict,
+    embed_dicom_and_nifti_metadata,
+    group_dicoms_into_seqinfos,
+    parse_private_csa_header,
+)
 from .utils import (
     assert_cwd_unchanged,
     TESTS_DATA_PATH,
@@ -64,3 +70,18 @@ def test_embed_dicom_and_nifti_metadata(tmpdir):
 
     assert out3.pop("existing") == "data"
     assert out3 == out2
+
+
+def test_group_dicoms_into_seqinfos(tmpdir):
+    """Tests for group_dicoms_into_seqinfos"""
+
+    # 1) Check that it works for PhoenixDocuments:
+    # set up testing files
+    dcmfolder = op.join(TESTS_DATA_PATH, 'Phoenix')
+    dcmfiles = glob(op.join(dcmfolder, '*', '*.dcm'))
+
+    seqinfo = group_dicoms_into_seqinfos(dcmfiles, 'studyUID', flatten=True)
+
+    assert type(seqinfo) is OrderedDict
+    assert len(seqinfo) == len(dcmfiles)
+    assert [s.series_description for s in seqinfo] == ['AAHead_Scout_32ch-head-coil', 'PhoenixZIPReport']

heudiconv/tests/test_heuristics.py

@@ -176,3 +176,20 @@ def test_notop(tmpdir, bidsoptions):
             assert not op.exists(pjoin(tmppath, 'Halchenko/Yarik/950_bids_test4', fname))
         else:
             assert op.exists(pjoin(tmppath, 'Halchenko/Yarik/950_bids_test4', fname))
+
+
+def test_phoenix_doc_conversion(tmpdir):
+    tmppath = tmpdir.strpath
+    subID = 'Phoenix'
+    args = (
+        "-c dcm2niix -o %s -b -f bids_PhoenixReport --files %s -s %s"
+        % (tmpdir, pjoin(TESTS_DATA_PATH, 'Phoenix'), subID)
+    ).split(' ')
+    runner(args)
+
+    # check that the Phoenix document has been extracted (as gzipped dicom) in
+    # the sourcedata/misc folder:
+    assert op.exists(pjoin(tmppath, 'sourcedata', 'sub-%s', 'misc', 'sub-%s_phoenix.dicom.tgz') % (subID, subID))
+    # check that no "sub-<subID>/misc" folder has been created in the BIDS
+    # structure:
+    assert not op.exists(pjoin(tmppath, 'sub-%s', 'misc') % subID)

heudiconv/tests/test_main.py

@@ -283,6 +283,11 @@ def test_cache(tmpdir):
     assert (cachedir / 'S01.auto.txt').exists()
     assert (cachedir / 'S01.edit.txt').exists()
 
+    # check dicominfo has "time" as last column:
+    with open(str(cachedir / 'dicominfo.tsv'), 'r') as f:
+        cols = f.readline().split()
+    assert cols[26] == "time"
+
 
 def test_no_etelemetry():
     # smoke test at large - just verifying that no crash if no etelemetry

heudiconv/tests/test_regression.py

@@ -9,14 +9,25 @@
 from heudiconv.external.pydicom import dcm
 from heudiconv.utils import load_json
 # testing utilities
-from .utils import fetch_data, gen_heudiconv_args, TESTS_DATA_PATH
+from .utils import (
+    fetch_data,
+    gen_heudiconv_args,
+    assert_cwd_unchanged,
+    TESTS_DATA_PATH
+)
 
 have_datalad = True
 try:
     from datalad.support.exceptions import IncompleteResultsError
 except ImportError:
     have_datalad = False
 
+have_bidsphysio = True
+try:
+    from bidsphysio.dcm2bids.dcm2bidsphysio import dcm2bids
+except ImportError:
+    have_bidsphysio = False
+
 
 @pytest.mark.skipif(not have_datalad, reason="no datalad")
 @pytest.mark.parametrize('subject', ['sub-sid000143'])
@@ -84,6 +95,33 @@ def test_multiecho(tmpdir, subject='MEEPI', heuristic='bids_ME.py'):
         assert 'echo-' not in event
 
 
+@assert_cwd_unchanged(ok_to_chdir=True)  # so we cd back after tmpdir.chdir
+@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 = "{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 = [

heudiconv/utils.py

@@ -46,7 +46,8 @@
     'patient_sex',           # 23
     'date',                  # 24
     'series_uid',            # 25
- ]
+    'time',                  # 26
+]
 
 SeqInfo = namedtuple('SeqInfo', seqinfo_fields)