[ENH] Add phase1/phase2 SDC support

.zenodo.json

@@ -84,6 +84,10 @@
       "orcid": "0000-0002-8508-9088"
     },
     {
+      "name": "Cieslak, Matthew",
+      "affiliation": "Department of Neuropsychiatry, University of Pennsylvania",
+      "orcid": "0000-0002-1931-4734"
+    },
       "name": "Liem, Franz",
       "affiliation": "University of Zurich",
       "orcid": "0000-0003-0646-4810"

fmriprep/interfaces/__init__.py

@@ -2,7 +2,8 @@
 # emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
 # vi: set ft=python sts=4 ts=4 sw=4 et:
 from .bids import (
-    ReadSidecarJSON, DerivativesDataSink, BIDSDataGrabber, BIDSFreeSurferDir, BIDSInfo
+    ReadSidecarJSON, DerivativesDataSink, BIDSDataGrabber, BIDSFreeSurferDir,
+    BIDSInfo
 )
 from .images import (
     IntraModalMerge, ValidateImage, TemplateDimensions, Conform, MatchHeader
@@ -13,8 +14,10 @@
 )
 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 .utils import (TPM2ROI, AddTPMs, AddTSVHeader, ConcatAffines,
+                    JoinTSVColumns)
+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

@@ -19,9 +19,8 @@
 import nibabel as nb
 from nipype import logging
 from nipype.utils.filemanip import fname_presuffix
-from nipype.interfaces.base import (
-    BaseInterfaceInputSpec, TraitedSpec, File, isdefined, traits,
-    SimpleInterface)
+from nipype.interfaces.base import (BaseInterfaceInputSpec, TraitedSpec, File, isdefined, traits,
+                                    SimpleInterface, InputMultiObject)
 
 LOGGER = logging.getLogger('nipype.interface')
 
@@ -69,10 +68,9 @@ def _run_interface(self, runtime):
             # Dilate mask
             if self.inputs.mask_erode > 0:
                 struc = sim.iterate_structure(sim.generate_binary_structure(3, 2), 1)
+                # pylint: disable=no-member
                 mask = sim.binary_erosion(
-                    mask, struc,
-                    iterations=self.inputs.mask_erode
-                ).astype(np.uint8)  # pylint: disable=no-member
+                    mask, struc, iterations=self.inputs.mask_erode).astype(np.uint8)
 
         self._results['out_file'] = fname_presuffix(
             self.inputs.in_file, suffix='_enh', newpath=runtime.cwd)
@@ -82,8 +80,8 @@ def _run_interface(self, runtime):
             data = _unwrap(data, self.inputs.in_magnitude, mask)
             self._results['out_unwrapped'] = fname_presuffix(
                 self.inputs.in_file, suffix='_unwrap', newpath=runtime.cwd)
-            nb.Nifti1Image(data, fmap_nii.affine, fmap_nii.header).to_filename(
-                self._results['out_unwrapped'])
+            nb.Nifti1Image(data, fmap_nii.affine,
+                           fmap_nii.header).to_filename(self._results['out_unwrapped'])
 
         if not self.inputs.bspline_smooth:
             datanii.to_filename(self._results['out_file'])
@@ -93,8 +91,7 @@ def _run_interface(self, runtime):
             from statsmodels.robust.scale import mad
 
             # Fit BSplines (coarse)
-            bspobj = fbsp.BSplineFieldmap(datanii, weights=mask,
-                                          njobs=self.inputs.num_threads)
+            bspobj = fbsp.BSplineFieldmap(datanii, weights=mask, njobs=self.inputs.num_threads)
             bspobj.fit()
             smoothed1 = bspobj.get_smoothed()
 
@@ -116,11 +113,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 +126,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,9 +198,34 @@ class Phasediff2Fieldmap(SimpleInterface):
 
     def _run_interface(self, runtime):
         self._results['out_file'] = phdiff2fmap(
-            self.inputs.in_file,
-            _delta_te(self.inputs.metadata),
-            newpath=runtime.cwd)
+            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
 
 
@@ -233,8 +255,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
 
@@ -255,9 +276,10 @@ def _unwrap(fmap_data, mag_file, mask=None):
     nb.Nifti1Image(magnii.get_data(), magnii.affine).to_filename('fmap_mag.nii.gz')
 
     # Run prelude
-    res = PRELUDE(phase_file='fmap_rad.nii.gz',
-                  magnitude_file='fmap_mag.nii.gz',
-                  mask_file='fmap_mask.nii.gz').run()
+    res = PRELUDE(
+        phase_file='fmap_rad.nii.gz',
+        magnitude_file='fmap_mag.nii.gz',
+        mask_file='fmap_mask.nii.gz').run()
 
     unwrapped = nb.load(res.outputs.unwrapped_phase_file).get_data() * (fmapmax / pi)
     return unwrapped
@@ -290,8 +312,9 @@ def get_ees(in_meta, in_file=None):
 
            =  T_\\text{ro} \\,  (N_\\text{PE} / f_\\text{acc} - 1)^{-1}
 
-    where :math:`N_y` is the number of pixels along the phase-encoding direction
-    :math:`y`, and :math:`f_\\text{acc}` is the parallel imaging acceleration factor
+    where :math:`N_y` is the number of pixels along the phase-encoding
+    direction :math:`y`, and :math:`f_\\text{acc}` is the parallel imaging
+    acceleration factor
     (:abbr:`GRAPPA (GeneRalized Autocalibrating Partial Parallel Acquisition)`,
     :abbr:`ARC (Autocalibrating Reconstruction for Cartesian imaging)`, etc.).
 
@@ -471,7 +494,8 @@ 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 +521,51 @@ 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()
+
+    # Calculate fieldmaps
+    rad0 = np.pi * phase0 / 2048 - np.pi
+    rad1 = np.pi * phase1 / 2048 - np.pi
+    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):
@@ -520,13 +589,13 @@ def _delta_te(in_values, te1=None, te2=None):
 
     # For convienience if both are missing we should give one error about them
     if te1 is None and te2 is None:
-        raise RuntimeError('EchoTime1 and EchoTime2 metadata fields not found. '
-                           'Please consult the BIDS specification.')
+        raise RuntimeError('EchoTime1 and EchoTime2 metadata fields not found.'
+                           ' Please consult the BIDS specification.')
     if te1 is None:
-        raise RuntimeError(
-            'EchoTime1 metadata field not found. Please consult the BIDS specification.')
+        raise RuntimeError('EchoTime1 metadata field not found. Please consult the BIDS '
+                           'specification.')
     if te2 is None:
-        raise RuntimeError(
-            'EchoTime2 metadata field not found. Please consult the BIDS specification.')
+        raise RuntimeError('EchoTime2 metadata field not found. Please consult the BIDS '
+                           'specification.')
 
     return abs(float(te2) - float(te1))

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

@@ -9,8 +9,8 @@
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 If the dataset metadata indicate tha more than one field map acquisition is
-``IntendedFor`` (see BIDS Specification section 8.9) the following priority will
-be used:
+``IntendedFor`` (see BIDS Specification section 8.9) the following priority
+will be used:
 
   1. :ref:`sdc_pepolar` (or **blip-up/blip-down**)
 
@@ -52,7 +52,8 @@
     'epi': 0,
     'fieldmap': 1,
     'phasediff': 2,
-    'syn': 3
+    'phase': 3,
+    'syn': 4
 }
 DEFAULT_MEMORY_MIN_GB = 0.01
 
@@ -80,13 +81,17 @@ def init_sdc_wf(fmaps, bold_meta, omp_nthreads=1,
         wf = init_sdc_wf(
             fmaps=[{
                 'type': 'phasediff',
-                'phasediff': 'sub-03/ses-2/fmap/sub-03_ses-2_run-1_phasediff.nii.gz',
-                'magnitude1': 'sub-03/ses-2/fmap/sub-03_ses-2_run-1_magnitude1.nii.gz',
-                'magnitude2': 'sub-03/ses-2/fmap/sub-03_ses-2_run-1_magnitude2.nii.gz',
+                'phasediff': \
+                    'sub-03/ses-2/fmap/sub-03_ses-2_run-1_phasediff.nii.gz',
+                'magnitude1': \
+                    'sub-03/ses-2/fmap/sub-03_ses-2_run-1_magnitude1.nii.gz',
+                'magnitude2': \
+                    'sub-03/ses-2/fmap/sub-03_ses-2_run-1_magnitude2.nii.gz',
             }],
             bold_meta={
                 'RepetitionTime': 2.0,
-                'SliceTiming': [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9],
+                'SliceTiming': [0.0, 0.1, 0.2, 0.3, 0.4, 0.5,
+                                0.6, 0.7, 0.8, 0.9],
                 'PhaseEncodingDirection': 'j',
             },
         )
@@ -175,10 +180,12 @@ def init_sdc_wf(fmaps, bold_meta, omp_nthreads=1,
 
     # PEPOLAR path
     if fmap['type'] == 'epi':
-        outputnode.inputs.method = 'PEB/PEPOLAR (phase-encoding based / PE-POLARity)'
+        outputnode.inputs.method = \
+            'PEB/PEPOLAR (phase-encoding based / PE-POLARity)'
         # Get EPI polarities and their metadata
-        epi_fmaps = [(fmap_['epi'], fmap_['metadata']["PhaseEncodingDirection"])
-                     for fmap_ in fmaps if fmap_['type'] == 'epi']
+        epi_fmaps = [
+            (fmap_['epi'], fmap_['metadata']["PhaseEncodingDirection"])
+            for fmap_ in fmaps if fmap_['type'] == 'epi']
         sdc_unwarp_wf = init_pepolar_unwarp_wf(
             bold_meta=bold_meta,
             epi_fmaps=epi_fmaps,
@@ -193,7 +200,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 +213,18 @@ 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")
@@ -243,7 +257,8 @@ def init_sdc_wf(fmaps, bold_meta, omp_nthreads=1,
         workflow.connect([
             (inputnode, syn_sdc_wf, [
                 ('t1_brain', 'inputnode.t1_brain'),
-                ('t1_2_mni_reverse_transform', 'inputnode.t1_2_mni_reverse_transform'),
+                ('t1_2_mni_reverse_transform',
+                    'inputnode.t1_2_mni_reverse_transform'),
                 ('bold_ref', 'inputnode.bold_ref'),
                 ('bold_ref_brain', 'inputnode.bold_ref_brain'),
                 ('template', 'inputnode.template')]),