[ENH] Add phase1/phase2 SDC support

.zenodo.json

@@ -106,6 +106,11 @@
       "affiliation": "University College London",
       "orcid": "0000-0002-9699-9052"
     },
+    {
+      "name": "Cieslak, Matthew",
+      "affiliation": "Department of Neuropsychiatry, University of Pennsylvania",
+      "orcid": "0000-0002-1931-4734"
+    },
     {
       "name": "Halchenko, Yaroslav O.",
       "affiliation": "Dartmouth College: Hanover, NH, United States",

fmriprep/data/boilerplate.bib

@@ -306,6 +306,16 @@ @article{hcppipelines
     year = 2013
 }
 
+@article{pncprocessing,
+  title={Neuroimaging of the Philadelphia neurodevelopmental cohort},
+  author={Satterthwaite, Theodore D and Elliott, Mark A and Ruparel, Kosha and Loughead, James and Prabhakaran, Karthik and Calkins, Monica E and Hopson, Ryan and Jackson, Chad and Keefe, Jack and Riley, Marisa and others},
+  journal={Neuroimage},
+  volume={86},
+  pages={544--553},
+  year={2014},
+  publisher={Elsevier}
+}
+
 @article{fs_template,
     author = {Reuter, Martin and Rosas, Herminia Diana and Fischl, Bruce},
     doi = {10.1016/j.neuroimage.2010.07.020},

fmriprep/interfaces/__init__.py

@@ -6,7 +6,7 @@
     bids, cifti, freesurfer, images, itk, surf, utils)
 
 from .reports import SubjectSummary, FunctionalSummary, AboutSummary
-from .fmap import FieldEnhance, FieldToRadS, FieldToHz, Phasediff2Fieldmap
+from .fmap import FieldEnhance, FieldToRadS, FieldToHz, Phasediff2Fieldmap, Phases2Fieldmap
 from .confounds import GatherConfounds, ICAConfounds, FMRISummary
 from .multiecho import T2SMap
 

fmriprep/interfaces/fmap.py

@@ -21,7 +21,7 @@
 from nipype.utils.filemanip import fname_presuffix
 from nipype.interfaces.base import (
     BaseInterfaceInputSpec, TraitedSpec, File, isdefined, traits,
-    SimpleInterface)
+    SimpleInterface, InputMultiObject)
 
 LOGGER = logging.getLogger('nipype.interface')
 
@@ -207,6 +207,34 @@ def _run_interface(self, runtime):
         return runtime
 
 
+class Phases2FieldmapInputSpec(BaseInterfaceInputSpec):
+    phase_files = InputMultiObject(
+        File(exists=True), mandatory=True, desc='list of phase1, phase2 files')
+    metadatas = traits.List(
+        traits.Dict, mandatory=True, desc='list of phase1, phase2 metadata dicts')
+
+
+class Phases2FieldmapOutputSpec(TraitedSpec):
+    out_file = File(desc='the output fieldmap')
+    phasediff_metadata = traits.Dict(desc='the phasediff metadata')
+
+
+class Phases2Fieldmap(SimpleInterface):
+    """
+    Convert a phase1, phase2 into a difference map
+    """
+    input_spec = Phases2FieldmapInputSpec
+    output_spec = Phases2FieldmapOutputSpec
+
+    def _run_interface(self, runtime):
+        # Get the echo times
+        fmap_file, merged_metadata = phases2fmap(self.inputs.phase_files, self.inputs.metadatas,
+                                                 newpath=runtime.cwd)
+        self._results['phasediff_metadata'] = merged_metadata
+        self._results['out_file'] = fmap_file
+        return runtime
+
+
 def _despike2d(data, thres, neigh=None):
     """
     despiking as done in FSL fugue
@@ -497,8 +525,69 @@ def phdiff2fmap(in_file, delta_te, newpath=None):
     return out_file
 
 
+def phases2fmap(phase_files, metadatas, newpath=None):
+    """Calculates a phasediff from two phase images. Assumes monopolar
+    readout. """
+    import numpy as np
+    import nibabel as nb
+    from nipype.utils.filemanip import fname_presuffix
+    from copy import deepcopy
+
+    phasediff_file = fname_presuffix(phase_files[0], suffix='_phasediff', newpath=newpath)
+    echo_times = [meta.get("EchoTime") for meta in metadatas]
+    if None in echo_times or echo_times[0] == echo_times[1]:
+        raise RuntimeError()
+    # Determine the order of subtraction
+    short_echo_index = echo_times.index(min(echo_times))
+    long_echo_index = echo_times.index(max(echo_times))
+
+    short_phase_image = phase_files[short_echo_index]
+    long_phase_image = phase_files[long_echo_index]
+
+    image0 = nb.load(short_phase_image)
+    phase0 = image0.get_fdata()
+    image1 = nb.load(long_phase_image)
+    phase1 = image1.get_fdata()
+
+    def rescale_image(img):
+        if np.any(img < -128):
+            # This happens sometimes on 7T fieldmaps
+            LOGGER.info("Found negative values in phase image: rescaling")
+            imax = img.max()
+            imin = img.min()
+            scaled = 2 * ((img - imin) / (imax - imin) - 0.5)
+            return np.pi * scaled
+        mask = img > 0
+        imax = img.max()
+        imin = img.min()
+        max_check = imax - 4096
+        if np.abs(max_check) > 10 or np.abs(imin) > 10:
+            LOGGER.warning("Phase image may be scaled incorrectly: check results")
+        return mask * (img / 2048 * np.pi - np.pi)
+
+    # Calculate fieldmaps
+    rad0 = rescale_image(phase0)
+    rad1 = rescale_image(phase1)
+    a = np.cos(rad0)
+    b = np.sin(rad0)
+    c = np.cos(rad1)
+    d = np.sin(rad1)
+    fmap = -np.arctan2(b * c - a * d, a * c + b * d)
+
+    phasediff_nii = nb.Nifti1Image(fmap, image0.affine)
+    phasediff_nii.set_data_dtype(np.float32)
+    phasediff_nii.to_filename(phasediff_file)
+
+    merged_metadata = deepcopy(metadatas[0])
+    del merged_metadata['EchoTime']
+    merged_metadata['EchoTime1'] = float(echo_times[short_echo_index])
+    merged_metadata['EchoTime2'] = float(echo_times[long_echo_index])
+
+    return phasediff_file, merged_metadata
+
+
 def _delta_te(in_values, te1=None, te2=None):
-    """Read :math:`\Delta_\text{TE}` from BIDS metadata dict"""
+    r"""Read :math:`\Delta_\text{TE}` from BIDS metadata dict"""
     if isinstance(in_values, float):
         te2 = in_values
         te1 = 0.

fmriprep/workflows/bold/base.py

@@ -336,7 +336,11 @@ def init_func_preproc_wf(
         if 'fieldmaps' not in ignore:
             fmaps = layout.get_fieldmap(ref_file, return_list=True)
             for fmap in fmaps:
-                fmap['metadata'] = layout.get_metadata(fmap[fmap['suffix']])
+                if fmap['suffix'] == 'phase':
+                    fmap_key = 'phase1'
+                else:
+                    fmap_key = fmap['suffix']
+                fmap['metadata'] = layout.get_metadata(fmap[fmap_key])
 
         # Run SyN if forced or in the absence of fieldmap correction
         if force_syn or (use_syn and not fmaps):

fmriprep/workflows/fieldmap/base.py

@@ -53,7 +53,8 @@
     'epi': 0,
     'fieldmap': 1,
     'phasediff': 2,
-    'syn': 3
+    'phase': 3,
+    'syn': 4
 }
 DEFAULT_MEMORY_MIN_GB = 0.01
 
@@ -199,7 +200,7 @@ def init_sdc_wf(fmaps, bold_meta, omp_nthreads=1,
         ])
 
     # FIELDMAP path
-    if fmap['suffix'] in ['fieldmap', 'phasediff']:
+    if fmap['suffix'] in ['fieldmap', 'phasediff', 'phase']:
         outputnode.inputs.method = 'FMB (%s-based)' % fmap['suffix']
         # Import specific workflows here, so we don't break everything with one
         # unused workflow.
@@ -212,11 +213,30 @@ def init_sdc_wf(fmaps, bold_meta, omp_nthreads=1,
             fmap_estimator_wf.inputs.inputnode.fieldmap = fmap['fieldmap']
             fmap_estimator_wf.inputs.inputnode.magnitude = fmap['magnitude']
 
-        if fmap['suffix'] == 'phasediff':
+        if fmap['suffix'] in ('phasediff', 'phase'):
             from .phdiff import init_phdiff_wf
-            fmap_estimator_wf = init_phdiff_wf(omp_nthreads=omp_nthreads)
+            fmap_estimator_wf = init_phdiff_wf(omp_nthreads=omp_nthreads,
+                                               phasetype=fmap['suffix'])
             # set inputs
-            fmap_estimator_wf.inputs.inputnode.phasediff = fmap['phasediff']
+            if fmap['suffix'] == 'phasediff':
+                fmap_estimator_wf.inputs.inputnode.phasediff = fmap['phasediff']
+            elif fmap['suffix'] == 'phase':
+                # Check that fieldmap is not bipolar
+                fmap_polarity = fmap['metadata'].get('DiffusionScheme', None)
+                if fmap_polarity == 'Bipolar':
+                    LOGGER.warning("Bipolar fieldmaps are not supported. Ignoring")
+                    workflow.__postdesc__ = ""
+                    outputnode.inputs.method = 'None'
+                    workflow.connect([
+                        (inputnode, outputnode, [('bold_ref', 'bold_ref'),
+                                                 ('bold_mask', 'bold_mask'),
+                                                 ('bold_ref_brain', 'bold_ref_brain')]),
+                    ])
+                    return workflow
+                if fmap_polarity is None:
+                    LOGGER.warning("Assuming phase images are Monopolar")
+
+                fmap_estimator_wf.inputs.inputnode.phasediff = [fmap['phase1'], fmap['phase2']]
             fmap_estimator_wf.inputs.inputnode.magnitude = [
                 fmap_ for key, fmap_ in sorted(fmap.items())
                 if key.startswith("magnitude")

fmriprep/workflows/fieldmap/phdiff.py

@@ -15,6 +15,7 @@
 
 Fieldmap preprocessing workflow for fieldmap data structure
 8.9.1 in BIDS 1.0.0: one phase diff and at least one magnitude image
+8.9.2 in BIDS 1.0.0: two phases and at least one magnitude image
 
 """
 
@@ -27,10 +28,10 @@
 from niworkflows.interfaces.images import IntraModalMerge
 from niworkflows.interfaces.masks import BETRPT
 
-from ...interfaces import Phasediff2Fieldmap, DerivativesDataSink
+from ...interfaces import Phasediff2Fieldmap, Phases2Fieldmap, DerivativesDataSink
 
 
-def init_phdiff_wf(omp_nthreads, name='phdiff_wf'):
+def init_phdiff_wf(omp_nthreads, phasetype='phasediff', name='phdiff_wf'):
     """
     Estimates the fieldmap using a phase-difference image and one or more
     magnitude images corresponding to two or more :abbr:`GRE (Gradient Echo sequence)`
@@ -69,12 +70,6 @@ def init_phdiff_wf(omp_nthreads, name='phdiff_wf'):
     outputnode = pe.Node(niu.IdentityInterface(
         fields=['fmap', 'fmap_ref', 'fmap_mask']), name='outputnode')
 
-    def _pick1st(inlist):
-        return inlist[0]
-
-    # Read phasediff echo times
-    meta = pe.Node(ReadSidecarJSON(bids_validate=False), name='meta', mem_gb=0.01)
-
     # Merge input magnitude images
     magmrg = pe.Node(IntraModalMerge(), name='magmrg')
 
@@ -90,9 +85,6 @@ def _pick1st(inlist):
     # dilate = pe.Node(fsl.maths.MathsCommand(
     #     nan2zeros=True, args='-kernel sphere 5 -dilM'), name='MskDilate')
 
-    # phase diff -> radians
-    pha2rads = pe.Node(niu.Function(function=siemens2rads), name='pha2rads')
-
     # FSL PRELUDE will perform phase-unwrapping
     prelude = pe.Node(fsl.PRELUDE(), name='prelude')
 
@@ -110,16 +102,47 @@ def _pick1st(inlist):
     # pre_fugue = pe.Node(fsl.FUGUE(save_fmap=True), name='ComputeFieldmapFUGUE')
     # rsec2hz (divide by 2pi)
 
+    if phasetype == "phasediff":
+        # Read phasediff echo times
+        meta = pe.Node(ReadSidecarJSON(bids_validate=False), name='meta', mem_gb=0.01)
+
+        # phase diff -> radians
+        pha2rads = pe.Node(niu.Function(function=siemens2rads),
+                           name='pha2rads')
+        # Read phasediff echo times
+        meta = pe.Node(ReadSidecarJSON(), name='meta', mem_gb=0.01,
+                       run_without_submitting=True)
+        workflow.connect([
+            (meta, compfmap, [('out_dict', 'metadata')]),
+            (inputnode, pha2rads, [('phasediff', 'in_file')]),
+            (pha2rads, prelude, [('out', 'phase_file')]),
+            (inputnode, ds_report_fmap_mask, [('phasediff', 'source_file')]),
+        ])
+
+    elif phasetype == "phase":
+        workflow.__desc__ += """\
+The phase difference used for unwarping was calculated using two separate phase measurements
+ [@pncprocessing].
+    """
+        # Special case for phase1, phase2 images
+        meta = pe.MapNode(ReadSidecarJSON(), name='meta', mem_gb=0.01,
+                          run_without_submitting=True, iterfield=['in_file'])
+        phases2fmap = pe.Node(Phases2Fieldmap(), name='phases2fmap')
+        workflow.connect([
+            (meta, phases2fmap, [('out_dict', 'metadatas')]),
+            (inputnode, phases2fmap, [('phasediff', 'phase_files')]),
+            (phases2fmap, prelude, [('out_file', 'phase_file')]),
+            (phases2fmap, compfmap, [('phasediff_metadata', 'metadata')]),
+            (phases2fmap, ds_report_fmap_mask, [('out_file', 'source_file')])
+        ])
+
     workflow.connect([
         (inputnode, meta, [('phasediff', 'in_file')]),
         (inputnode, magmrg, [('magnitude', 'in_files')]),
         (magmrg, n4, [('out_avg', 'input_image')]),
         (n4, prelude, [('output_image', 'magnitude_file')]),
         (n4, bet, [('output_image', 'in_file')]),
         (bet, prelude, [('mask_file', 'mask_file')]),
-        (inputnode, pha2rads, [('phasediff', 'in_file')]),
-        (pha2rads, prelude, [('out', 'phase_file')]),
-        (meta, compfmap, [('out_dict', 'metadata')]),
         (prelude, denoise, [('unwrapped_phase_file', 'in_file')]),
         (denoise, demean, [('out_file', 'in_file')]),
         (demean, cleanup_wf, [('out', 'inputnode.in_file')]),
@@ -128,7 +151,6 @@ def _pick1st(inlist):
         (compfmap, outputnode, [('out_file', 'fmap')]),
         (bet, outputnode, [('mask_file', 'fmap_mask'),
                            ('out_file', 'fmap_ref')]),
-        (inputnode, ds_report_fmap_mask, [('phasediff', 'source_file')]),
         (bet, ds_report_fmap_mask, [('out_report', 'in_file')]),
     ])