[ENH] Add phase1/phase2 SDC support

fmriprep/interfaces/fmap.py

@@ -207,6 +207,27 @@ def _run_interface(self, runtime):
         return runtime
 
 
+class Phases2FieldmapInputSpec(BaseInterfaceInputSpec):
+    in_file = File(exists=True, mandatory=True, desc='input fieldmap')
+    metadata = traits.List(mandatory=True, desc='BIDS metadata dictionaries')
+
+
+class Phases2Fieldmap(SimpleInterface):
+    """
+    Convert a phase difference map into a fieldmap in Hz
+    """
+    input_spec = Phases2FieldmapInputSpec
+    output_spec = Phasediff2FieldmapOutputSpec
+
+    def _run_interface(self, runtime):
+        self._results['out_file'] = phdiff2fmap(
+            self.inputs.in_file,
+            _delta_te_from_two_phases(self.inputs.metadata),
+            newpath=runtime.cwd)
+        return runtime
+
+
+
 def _despike2d(data, thres, neigh=None):
     """
     despiking as done in FSL fugue
@@ -496,6 +517,39 @@ def phdiff2fmap(in_file, delta_te, newpath=None):
     nii.to_filename(out_file)
     return out_file
 
+def phases2fmap(in_files, delta_te, newpath=None):
+    r"""
+    Converts the input phase maps (with different echo times) into a fieldmap
+    in Hz, using the eq. (1) of [Hutton2002]_:
+
+    .. math::
+
+        \Delta B_0 (\text{T}^{-1}) = \frac{\Delta \Theta}{2\pi\gamma \Delta\text{TE}}
+
+
+    In this case, we do not take into account the gyromagnetic ratio of the
+    proton (:math:`\gamma`), since it will be applied inside TOPUP:
+
+    .. math::
+
+        \Delta B_0 (\text{Hz}) = \frac{\Delta \Theta}{2\pi \Delta\text{TE}}
+
+    """
+    import math
+    import numpy as np
+    import nibabel as nb
+    from nipype.utils.filemanip import fname_presuffix
+    #  GYROMAG_RATIO_H_PROTON_MHZ = 42.576
+
+
+    out_file = fname_presuffix(in_file, suffix='_fmap', newpath=newpath)
+    image = nb.load(in_file)
+    data = (image.get_data().astype(np.float32) / (2. * math.pi * delta_te))
+    nii = nb.Nifti1Image(data, image.affine, image.header)
+    nii.set_data_dtype(np.float32)
+    nii.to_filename(out_file)
+    return out_file
+
 
 def _delta_te(in_values, te1=None, te2=None):
     """Read :math:`\Delta_\text{TE}` from BIDS metadata dict"""
@@ -530,3 +584,30 @@ def _delta_te(in_values, te1=None, te2=None):
             'EchoTime2 metadata field not found. Please consult the BIDS specification.')
 
     return abs(float(te2) - float(te1))
+
+def _delta_te_from_two_phases(in_dicts, te1=None, te2=None):
+    """Read :math:`\Delta_\text{TE}` from two BIDS metadata dicts"""
+
+    if isinstance(in_values, list):
+        te2_value, te1_value = in_values
+        if isinstance(te1_value, list):
+            te1_dict = te1[1]
+            if isinstance(te1_dict, dict):
+                te1 = te1_dict.get('EchoTime')
+        if isinstance(te2_value, list):
+            te2_dict = te2[1]
+            if isinstance(te2_dict, dict):
+                te2 = te2_dict.get('EchoTime')
+
+    # For convienience if both are missing we should give one error about them
+    if te1 is None and te2 is None:
+        raise RuntimeError('EchoTime metadata fields not found. '
+                           'Please consult the BIDS specification.')
+    if te1 is None:
+        raise RuntimeError(
+            'EchoTime metadata field not found for phase1. Please consult the BIDS specification.')
+    if te2 is None:
+        raise RuntimeError(
+            'EchoTime metadata field not found for phase2. Please consult the BIDS specification.')
+
+    return abs(float(te2) - float(te1))

fmriprep/workflows/fieldmap/base.py

@@ -52,7 +52,8 @@
     'epi': 0,
     'fieldmap': 1,
     'phasediff': 2,
-    'syn': 3
+    'phase':3,
+    'syn': 4
 }
 DEFAULT_MEMORY_MIN_GB = 0.01
 
@@ -193,7 +194,7 @@ def init_sdc_wf(fmaps, bold_meta, omp_nthreads=1,
         ])
 
     # FIELDMAP path
-    if fmap['type'] in ['fieldmap', 'phasediff']:
+    if fmap['type'] in ['fieldmap', 'phasediff', 'phase']:
         outputnode.inputs.method = 'FMB (%s-based)' % fmap['type']
         # Import specific workflows here, so we don't break everything with one
         # unused workflow.
@@ -206,11 +207,17 @@ 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['type'] == 'phasediff':
+        if fmap['type'] 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['type'])
             # set inputs
-            fmap_estimator_wf.inputs.inputnode.phasediff = fmap['phasediff']
+            if fmap['type'] == 'phasediff':
+                fmap_estimator_wf.inputs.inputnode.phasediff = fmap['phasediff']
+            elif fmap['type'] == 'phase':
+                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

@@ -27,11 +27,11 @@
 from ...engine import Workflow
 from ...interfaces import (
     ReadSidecarJSON, IntraModalMerge, DerivativesDataSink,
-    Phasediff2Fieldmap
+    Phasediff2Fieldmap, Phases2Fieldmap
 )
 
 
-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)`
@@ -73,9 +73,19 @@ def init_phdiff_wf(omp_nthreads, name='phdiff_wf'):
     def _pick1st(inlist):
         return inlist[0]
 
-    # Read phasediff echo times
-    meta = pe.Node(ReadSidecarJSON(), name='meta', mem_gb=0.01, run_without_submitting=True)
+    if phasetype == "phasediff":
+        # Read phasediff echo times
+        meta = pe.Node(ReadSidecarJSON(), name='meta', mem_gb=0.01,
+                        run_without_submitting=True)
+        compfmap = pe.Node(Phasediff2Fieldmap(), name='compfmap')
 
+    elif phasetype == "phase":
+        meta = pe.MapNode(ReadSidecarJSON(), name='meta', mem_gb=0.01,
+                        run_without_submitting=True,iterfield=['in_file'])
+        compfmap = pe.Node(Phases2Fieldmap(), name='compfmap')
+
+    pha2rads = pe.Node(niu.Function(function=siemens2rads), name='pha2rads')
+
     # Merge input magnitude images
     magmrg = pe.Node(IntraModalMerge(), name='magmrg')
 
@@ -90,9 +100,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')
 
@@ -103,7 +110,6 @@ def _pick1st(inlist):
 
     cleanup_wf = cleanup_edge_pipeline(name="cleanup_wf")
 
-    compfmap = pe.Node(Phasediff2Fieldmap(), name='compfmap')
 
     # The phdiff2fmap interface is equivalent to:
     # rad2rsec (using rads2radsec from nipype.workflows.dmri.fsl.utils)
@@ -118,8 +124,8 @@ def _pick1st(inlist):
         (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')]),
+        (pha2rads, prelude, [('out', 'phase_file')]),
         (prelude, denoise, [('unwrapped_phase_file', 'in_file')]),
         (denoise, demean, [('out_file', 'in_file')]),
         (demean, cleanup_wf, [('out', 'inputnode.in_file')]),