[ENH] Add phase1/phase2 SDC support

.zenodo.json

@@ -98,6 +98,11 @@
       "affiliation": "Department of Psychology, Columbia University",
       "orcid": "0000-0001-7936-9991"
     },
+    {
+      "name": "Cieslak, Matthew",
+      "affiliation": "Department of Neuropsychiatry, University of Pennsylvania",
+      "orcid": "0000-0002-1931-4734"
+    },
     {
       "name": "Poldrack, Russell A.",
       "affiliation": "Department of Psychology, Stanford University",

fmriprep/data/boilerplate.bib

@@ -271,6 +271,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

@@ -14,7 +14,7 @@
 from .surf import NormalizeSurf, GiftiNameSource, GiftiSetAnatomicalStructure
 from .reports import SubjectSummary, FunctionalSummary, AboutSummary
 from .utils import TPM2ROI, AddTPMs, AddTSVHeader, ConcatAffines, JoinTSVColumns
-from .fmap import FieldEnhance, FieldToRadS, FieldToHz, Phasediff2Fieldmap
+from .fmap import FieldEnhance, FieldToRadS, FieldToHz, Phasediff2Fieldmap, Phases2Fieldmap
 from .confounds import GatherConfounds, ICAConfounds, FMRISummary
 from .itk import MCFLIRT2ITK, MultiApplyTransforms
 from .multiecho import FirstEcho

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')
 
@@ -116,11 +116,11 @@ def _run_interface(self, runtime):
             if nslices > 1:
                 diffmapmsk = mask[..., errorslice[0]:errorslice[-1]]
                 diffmapnii = nb.Nifti1Image(
-                    diffmap[..., errorslice[0]:errorslice[-1]] * diffmapmsk,
-                    datanii.affine, datanii.header)
+                    diffmap[..., errorslice[0]:errorslice[-1]] * diffmapmsk, datanii.affine,
+                    datanii.header)
 
-                bspobj2 = fbsp.BSplineFieldmap(diffmapnii, knots_zooms=[24., 24., 4.],
-                                               njobs=self.inputs.num_threads)
+                bspobj2 = fbsp.BSplineFieldmap(
+                    diffmapnii, knots_zooms=[24., 24., 4.], njobs=self.inputs.num_threads)
                 bspobj2.fit()
                 smoothed2 = bspobj2.get_smoothed().get_data()
 
@@ -129,8 +129,8 @@ def _run_interface(self, runtime):
             else:
                 final = smoothed1.get_data()
 
-            nb.Nifti1Image(final, datanii.affine, datanii.header).to_filename(
-                self._results['out_file'])
+            nb.Nifti1Image(final, datanii.affine,
+                           datanii.header).to_filename(self._results['out_file'])
 
         return runtime
 
@@ -201,12 +201,38 @@ class Phasediff2Fieldmap(SimpleInterface):
 
     def _run_interface(self, runtime):
         self._results['out_file'] = phdiff2fmap(
-            self.inputs.in_file,
-            _delta_te(self.inputs.metadata),
+            self.inputs.in_file, _delta_te(self.inputs.metadata),
             newpath=runtime.cwd)
         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)
+        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
@@ -233,8 +259,7 @@ def _despike2d(data, thres, neigh=None):
                 patch_range = vals.max() - vals.min()
                 patch_med = np.median(vals)
 
-                if (patch_range > 1e-6 and
-                        (abs(thisval - patch_med) / patch_range) > thres):
+                if (patch_range > 1e-6 and (abs(thisval - patch_med) / patch_range) > thres):
                     data[i, j, k] = patch_med
     return data
 
@@ -471,7 +496,7 @@ def phdiff2fmap(in_file, delta_te, newpath=None):
 
     .. math::
 
-        \Delta B_0 (\text{T}^{-1}) = \frac{\Delta \Theta}{2\pi\gamma \Delta\text{TE}}
+        \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
@@ -497,6 +522,57 @@ 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):
+        imax = img.max()
+        imin = img.min()
+        scaled = 2 * ((img - imin) / (imax - imin) - 0.5)
+        return np.pi * scaled
+
+    # 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"""
     if isinstance(in_values, float):

fmriprep/workflows/bold/base.py

@@ -290,7 +290,11 @@ def init_func_preproc_wf(bold_file, ignore, freesurfer,
         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['type']])
+                if fmap['type'] == 'phase':
+                    fmap_key = 'phase1'
+                else:
+                    fmap_key = fmap['type']
+                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

@@ -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,15 +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")
-            ]
+                fmap_ for key, fmap_ in sorted(fmap.items()) if key.startswith("magnitude")]
 
         sdc_unwarp_wf = init_sdc_unwarp_wf(
             omp_nthreads=omp_nthreads,

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 phase images and two magnitude images
 
 """
 
@@ -27,11 +28,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)`
@@ -70,12 +71,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(), name='meta', mem_gb=0.01, run_without_submitting=True)
-
     # 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,44 @@ def _pick1st(inlist):
     # pre_fugue = pe.Node(fsl.FUGUE(save_fmap=True), name='ComputeFieldmapFUGUE')
     # rsec2hz (divide by 2pi)
 
+    if phasetype == "phasediff":
+        # 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_fmap_mask, [('phasediff', 'source_file')]),
+        ])
+
+    elif phasetype == "phase":
+        workflow.__desc__ += """\
+The phase difference used for unwarping was calculated using a method developed
+by Mark Elliott [@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_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 +148,6 @@ def _pick1st(inlist):
         (compfmap, outputnode, [('out_file', 'fmap')]),
         (bet, outputnode, [('mask_file', 'fmap_mask'),
                            ('out_file', 'fmap_ref')]),
-        (inputnode, ds_fmap_mask, [('phasediff', 'source_file')]),
         (bet, ds_fmap_mask, [('out_report', 'in_file')]),
     ])