From b2a937ab5e3451875b6dbf6bf6cd9654629c0375 Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Thu, 1 Nov 2018 10:07:06 -0400
Subject: [PATCH 01/34] first pass at 2phases SDC pipeline

---
 fmriprep/interfaces/fmap.py           | 48 +++++++++++++++++++++++++++
 fmriprep/workflows/fieldmap/base.py   | 17 +++++++---
 fmriprep/workflows/fieldmap/phdiff.py | 24 ++++++++------
 3 files changed, 74 insertions(+), 15 deletions(-)

diff --git a/fmriprep/interfaces/fmap.py b/fmriprep/interfaces/fmap.py
index f4286c3ea..6d633495a 100644
--- a/fmriprep/interfaces/fmap.py
+++ b/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
@@ -530,3 +551,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, list):
+            te1_dict = te1[1]
+            if isinstance(te1_dict, dict):
+                te1 = te1_dict.get('EchoTime')
+        if isinstance(te2, 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))
diff --git a/fmriprep/workflows/fieldmap/base.py b/fmriprep/workflows/fieldmap/base.py
index dfa974ddc..5135cca3c 100644
--- a/fmriprep/workflows/fieldmap/base.py
+++ b/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")
diff --git a/fmriprep/workflows/fieldmap/phdiff.py b/fmriprep/workflows/fieldmap/phdiff.py
index 23081889b..2e1a79e23 100644
--- a/fmriprep/workflows/fieldmap/phdiff.py
+++ b/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,8 +73,16 @@ 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')
 
     # Merge input magnitude images
     magmrg = pe.Node(IntraModalMerge(), name='magmrg')
@@ -90,9 +98,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 +108,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)
@@ -116,10 +120,10 @@ def _pick1st(inlist):
         (magmrg, n4, [('out_avg', 'input_image')]),
         (n4, prelude, [('output_image', 'magnitude_file')]),
         (n4, bet, [('output_image', 'in_file')]),
-        (bet, prelude, [('mask_file', 'mask_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')]),

From f2b45d4bc4baa36e98d36faee9cd6d2f669f9fdd Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Thu, 1 Nov 2018 10:26:31 -0400
Subject: [PATCH 02/34] workflow.connect syntax typo

---
 fmriprep/workflows/fieldmap/phdiff.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/fmriprep/workflows/fieldmap/phdiff.py b/fmriprep/workflows/fieldmap/phdiff.py
index 2e1a79e23..18d5d4da3 100644
--- a/fmriprep/workflows/fieldmap/phdiff.py
+++ b/fmriprep/workflows/fieldmap/phdiff.py
@@ -120,7 +120,7 @@ def _pick1st(inlist):
         (magmrg, n4, [('out_avg', 'input_image')]),
         (n4, prelude, [('output_image', 'magnitude_file')]),
         (n4, bet, [('output_image', 'in_file')]),
-        (bet, prelude, [('mask_file', 'mask_file')])]),
+        (bet, prelude, [('mask_file', 'mask_file')]),
         (inputnode, pha2rads, [('phasediff', 'in_file')]),
         (meta, compfmap, [('out_dict', 'metadata')]),
         (pha2rads, prelude, [('out', 'phase_file')]),

From 9dec5d8954414d584942b659d2e8507196fa807b Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Thu, 1 Nov 2018 10:45:34 -0400
Subject: [PATCH 03/34] remove unused variables from fmap.py

---
 fmriprep/interfaces/fmap.py | 37 +++++++++++++++++++++++++++++++++++--
 1 file changed, 35 insertions(+), 2 deletions(-)

diff --git a/fmriprep/interfaces/fmap.py b/fmriprep/interfaces/fmap.py
index 6d633495a..3e9fbe677 100644
--- a/fmriprep/interfaces/fmap.py
+++ b/fmriprep/interfaces/fmap.py
@@ -517,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"""
@@ -557,11 +590,11 @@ def _delta_te_from_two_phases(in_dicts, te1=None, te2=None):
 
     if isinstance(in_values, list):
         te2_value, te1_value = in_values
-        if isinstance(te1, list):
+        if isinstance(te1_value, list):
             te1_dict = te1[1]
             if isinstance(te1_dict, dict):
                 te1 = te1_dict.get('EchoTime')
-        if isinstance(te2, list):
+        if isinstance(te2_value, list):
             te2_dict = te2[1]
             if isinstance(te2_dict, dict):
                 te2 = te2_dict.get('EchoTime')

From 840d167869e3fc26d5d09cfe998954188612f11a Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Thu, 1 Nov 2018 10:50:58 -0400
Subject: [PATCH 04/34] add back pha2rads after accidental delete

---
 fmriprep/workflows/fieldmap/phdiff.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/fmriprep/workflows/fieldmap/phdiff.py b/fmriprep/workflows/fieldmap/phdiff.py
index 18d5d4da3..5f3bafea6 100644
--- a/fmriprep/workflows/fieldmap/phdiff.py
+++ b/fmriprep/workflows/fieldmap/phdiff.py
@@ -84,6 +84,8 @@ def _pick1st(inlist):
                         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')
 

From bb510de64ec09dd54cf4d6f3818cfac2b5ce3cd1 Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Thu, 1 Nov 2018 12:17:45 -0400
Subject: [PATCH 05/34] Changed copied code from phasediff functions

---
 fmriprep/interfaces/fmap.py           | 36 ++++++++++++++++-----------
 fmriprep/workflows/fieldmap/phdiff.py |  3 +--
 2 files changed, 22 insertions(+), 17 deletions(-)

diff --git a/fmriprep/interfaces/fmap.py b/fmriprep/interfaces/fmap.py
index 3e9fbe677..c1dc96d98 100644
--- a/fmriprep/interfaces/fmap.py
+++ b/fmriprep/interfaces/fmap.py
@@ -227,7 +227,6 @@ def _run_interface(self, runtime):
         return runtime
 
 
-
 def _despike2d(data, thres, neigh=None):
     """
     despiking as done in FSL fugue
@@ -517,6 +516,7 @@ 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
@@ -524,7 +524,7 @@ def phases2fmap(in_files, 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
@@ -541,11 +541,12 @@ def phases2fmap(in_files, delta_te, newpath=None):
     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)
+    out_file = fname_presuffix(in_files[0], suffix='_fmap', newpath=newpath)
+    image0 = nb.load(in_files[0])
+    image1 = nb.load(in_files[1])
+    data = image1.get_data() - image0.get_data()
+    data = (data / (2. * math.pi * delta_te))
+    nii = nb.Nifti1Image(data, image0.affine, image0.header)
     nii.set_data_dtype(np.float32)
     nii.to_filename(out_file)
     return out_file
@@ -574,22 +575,25 @@ 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.')
+            '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.')
+            '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(in_dicts, list):
+        te2_value, te1_value = in_dicts
         if isinstance(te1_value, list):
             te1_dict = te1[1]
             if isinstance(te1_dict, dict):
@@ -605,9 +609,11 @@ def _delta_te_from_two_phases(in_dicts, te1=None, te2=None):
                            'Please consult the BIDS specification.')
     if te1 is None:
         raise RuntimeError(
-            'EchoTime metadata field not found for phase1. Please consult the BIDS specification.')
+            '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.')
+            'EchoTime metadata field not found for phase2. Please consult the '
+            'BIDS specification.')
 
     return abs(float(te2) - float(te1))
diff --git a/fmriprep/workflows/fieldmap/phdiff.py b/fmriprep/workflows/fieldmap/phdiff.py
index 5f3bafea6..ec23bd435 100644
--- a/fmriprep/workflows/fieldmap/phdiff.py
+++ b/fmriprep/workflows/fieldmap/phdiff.py
@@ -85,7 +85,7 @@ def _pick1st(inlist):
         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')
 
@@ -110,7 +110,6 @@ def _pick1st(inlist):
 
     cleanup_wf = cleanup_edge_pipeline(name="cleanup_wf")
 
-
     # The phdiff2fmap interface is equivalent to:
     # rad2rsec (using rads2radsec from nipype.workflows.dmri.fsl.utils)
     # pre_fugue = pe.Node(fsl.FUGUE(save_fmap=True), name='ComputeFieldmapFUGUE')

From 4d7848334fceddfb9f66396444b61f39cf5e3d6a Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Thu, 1 Nov 2018 12:28:46 -0400
Subject: [PATCH 06/34] resolve flake8 errors

---
 fmriprep/interfaces/fmap.py         |  2 +-
 fmriprep/workflows/fieldmap/base.py | 34 ++++++++++++++++++-----------
 2 files changed, 22 insertions(+), 14 deletions(-)

diff --git a/fmriprep/interfaces/fmap.py b/fmriprep/interfaces/fmap.py
index c1dc96d98..cf48573c9 100644
--- a/fmriprep/interfaces/fmap.py
+++ b/fmriprep/interfaces/fmap.py
@@ -220,7 +220,7 @@ class Phases2Fieldmap(SimpleInterface):
     output_spec = Phasediff2FieldmapOutputSpec
 
     def _run_interface(self, runtime):
-        self._results['out_file'] = phdiff2fmap(
+        self._results['out_file'] = phases2fmap(
             self.inputs.in_file,
             _delta_te_from_two_phases(self.inputs.metadata),
             newpath=runtime.cwd)
diff --git a/fmriprep/workflows/fieldmap/base.py b/fmriprep/workflows/fieldmap/base.py
index 5135cca3c..e3d780925 100644
--- a/fmriprep/workflows/fieldmap/base.py
+++ b/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,7 @@
     'epi': 0,
     'fieldmap': 1,
     'phasediff': 2,
-    'phase':3,
+    'phase': 3,
     'syn': 4
 }
 DEFAULT_MEMORY_MIN_GB = 0.01
@@ -81,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',
             },
         )
@@ -176,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,
@@ -210,10 +216,11 @@ def init_sdc_wf(fmaps, bold_meta, omp_nthreads=1,
         if fmap['type'] in ('phasediff', 'phase'):
             from .phdiff import init_phdiff_wf
             fmap_estimator_wf = init_phdiff_wf(omp_nthreads=omp_nthreads,
-                                               phasetype = fmap['type'])
+                                               phasetype=fmap['type'])
             # set inputs
             if fmap['type'] == 'phasediff':
-                fmap_estimator_wf.inputs.inputnode.phasediff = fmap['phasediff']
+                fmap_estimator_wf.inputs.inputnode.phasediff = \
+                    fmap['phasediff']
             elif fmap['type'] == 'phase':
                 fmap_estimator_wf.inputs.inputnode.phasediff = [
                     fmap['phase1'], fmap['phase2']
@@ -250,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')]),

From 6e0eddaabe0a4ddfadcff2228e27630743fae9a5 Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Thu, 1 Nov 2018 13:14:32 -0400
Subject: [PATCH 07/34] More flake8 issues

---
 fmriprep/interfaces/fmap.py           | 34 +++++++++++++++++----------
 fmriprep/workflows/fieldmap/phdiff.py | 19 +++++++++------
 2 files changed, 34 insertions(+), 19 deletions(-)

diff --git a/fmriprep/interfaces/fmap.py b/fmriprep/interfaces/fmap.py
index cf48573c9..143e8ef26 100644
--- a/fmriprep/interfaces/fmap.py
+++ b/fmriprep/interfaces/fmap.py
@@ -32,10 +32,13 @@ class FieldEnhanceInputSpec(BaseInterfaceInputSpec):
     in_magnitude = File(exists=True, desc='input magnitude')
     unwrap = traits.Bool(False, usedefault=True, desc='run phase unwrap')
     despike = traits.Bool(True, usedefault=True, desc='run despike filter')
-    bspline_smooth = traits.Bool(True, usedefault=True, desc='run 3D bspline smoother')
+    bspline_smooth = traits.Bool(True, usedefault=True,
+                                 desc='run 3D bspline smoother')
     mask_erode = traits.Int(1, usedefault=True, desc='mask erosion iterations')
-    despike_threshold = traits.Float(0.2, usedefault=True, desc='mask erosion iterations')
-    num_threads = traits.Int(1, usedefault=True, nohash=True, desc='number of jobs')
+    despike_threshold = traits.Float(0.2, usedefault=True,
+                                     desc='mask erosion iterations')
+    num_threads = traits.Int(1, usedefault=True, nohash=True,
+                             desc='number of jobs')
 
 
 class FieldEnhanceOutputSpec(TraitedSpec):
@@ -68,7 +71,8 @@ def _run_interface(self, runtime):
 
             # Dilate mask
             if self.inputs.mask_erode > 0:
-                struc = sim.iterate_structure(sim.generate_binary_structure(3, 2), 1)
+                struc = sim.iterate_structure(
+                            sim.generate_binary_structure(3, 2), 1)
                 mask = sim.binary_erosion(
                     mask, struc,
                     iterations=self.inputs.mask_erode
@@ -108,7 +112,8 @@ def _run_interface(self, runtime):
 
             nslices = 0
             try:
-                errorslice = np.squeeze(np.argwhere(errormask.sum(0).sum(0) > 0))
+                errorslice = np.squeeze(
+                                np.argwhere(errormask.sum(0).sum(0) > 0))
                 nslices = errorslice[-1] - errorslice[0]
             except IndexError:  # mask is empty, do not refine
                 pass
@@ -119,7 +124,8 @@ def _run_interface(self, runtime):
                     diffmap[..., errorslice[0]:errorslice[-1]] * diffmapmsk,
                     datanii.affine, datanii.header)
 
-                bspobj2 = fbsp.BSplineFieldmap(diffmapnii, knots_zooms=[24., 24., 4.],
+                bspobj2 = fbsp.BSplineFieldmap(diffmapnii,
+                                               knots_zooms=[24., 24., 4.],
                                                njobs=self.inputs.num_threads)
                 bspobj2.fit()
                 smoothed2 = bspobj2.get_smoothed().get_data()
@@ -272,14 +278,16 @@ def _unwrap(fmap_data, mag_file, mask=None):
 
     nb.Nifti1Image(fmap_data, magnii.affine).to_filename('fmap_rad.nii.gz')
     nb.Nifti1Image(mask, magnii.affine).to_filename('fmap_mask.nii.gz')
-    nb.Nifti1Image(magnii.get_data(), magnii.affine).to_filename('fmap_mag.nii.gz')
+    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()
 
-    unwrapped = nb.load(res.outputs.unwrapped_phase_file).get_data() * (fmapmax / pi)
+    unwrapped = nb.load(res.outputs.unwrapped_phase_file
+                        ).get_data() * (fmapmax / pi)
     return unwrapped
 
 
@@ -310,8 +318,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.).
 
@@ -491,7 +500,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
@@ -616,4 +626,4 @@ def _delta_te_from_two_phases(in_dicts, te1=None, te2=None):
             'EchoTime metadata field not found for phase2. Please consult the '
             'BIDS specification.')
 
-    return abs(float(te2) - float(te1))
+    return float(te2) - float(te1)
diff --git a/fmriprep/workflows/fieldmap/phdiff.py b/fmriprep/workflows/fieldmap/phdiff.py
index ec23bd435..158ed3e8b 100644
--- a/fmriprep/workflows/fieldmap/phdiff.py
+++ b/fmriprep/workflows/fieldmap/phdiff.py
@@ -34,7 +34,8 @@
 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)`
+    magnitude images corresponding to two or more
+    :abbr:`GRE (Gradient Echo sequence)`
     acquisitions. The `original code was taken from nipype
     <https://github.com/nipy/nipype/blob/master/nipype/workflows/dmri/fsl/artifacts.py#L514>`_.
 
@@ -64,7 +65,8 @@ def init_phdiff_wf(omp_nthreads, phasetype="phasediff", name='phdiff_wf'):
 further improvements of HCP Pipelines [@hcppipelines].
 """
 
-    inputnode = pe.Node(niu.IdentityInterface(fields=['magnitude', 'phasediff']),
+    inputnode = pe.Node(niu.IdentityInterface(
+                        fields=['magnitude', 'phasediff']),
                         name='inputnode')
 
     outputnode = pe.Node(niu.IdentityInterface(
@@ -76,12 +78,12 @@ def _pick1st(inlist):
     if phasetype == "phasediff":
         # Read phasediff echo times
         meta = pe.Node(ReadSidecarJSON(), name='meta', mem_gb=0.01,
-                        run_without_submitting=True)
+                       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'])
+                          run_without_submitting=True, iterfield=['in_file'])
         compfmap = pe.Node(Phases2Fieldmap(), name='compfmap')
 
     pha2rads = pe.Node(niu.Function(function=siemens2rads), name='pha2rads')
@@ -94,7 +96,8 @@ def _pick1st(inlist):
                  name='n4', n_procs=omp_nthreads)
     bet = pe.Node(BETRPT(generate_report=True, frac=0.6, mask=True),
                   name='bet')
-    ds_fmap_mask = pe.Node(DerivativesDataSink(suffix='fmap_mask'), name='ds_report_fmap_mask',
+    ds_fmap_mask = pe.Node(DerivativesDataSink(suffix='fmap_mask'),
+                           name='ds_report_fmap_mask',
                            mem_gb=0.01, run_without_submitting=True)
     # uses mask from bet; outputs a mask
     # dilate = pe.Node(fsl.maths.MathsCommand(
@@ -103,7 +106,8 @@ def _pick1st(inlist):
     # FSL PRELUDE will perform phase-unwrapping
     prelude = pe.Node(fsl.PRELUDE(), name='prelude')
 
-    denoise = pe.Node(fsl.SpatialFilter(operation='median', kernel_shape='sphere',
+    denoise = pe.Node(fsl.SpatialFilter(operation='median',
+                                        kernel_shape='sphere',
                                         kernel_size=3), name='denoise')
 
     demean = pe.Node(niu.Function(function=demean_image), name='demean')
@@ -112,7 +116,8 @@ def _pick1st(inlist):
 
     # The phdiff2fmap interface is equivalent to:
     # rad2rsec (using rads2radsec from nipype.workflows.dmri.fsl.utils)
-    # pre_fugue = pe.Node(fsl.FUGUE(save_fmap=True), name='ComputeFieldmapFUGUE')
+    # pre_fugue = pe.Node(fsl.FUGUE(save_fmap=True),
+    #                     name='ComputeFieldmapFUGUE')
     # rsec2hz (divide by 2pi)
 
     workflow.connect([

From 9a0fcde8fe580356eccd8d3ffbb8a8ff9c1d5955 Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Thu, 1 Nov 2018 15:51:53 -0400
Subject: [PATCH 08/34] added import to interfaces/__init__.py

---
 fmriprep/interfaces/__init__.py | 2 +-
 fmriprep/workflows/bold/base.py | 6 +++++-
 2 files changed, 6 insertions(+), 2 deletions(-)

diff --git a/fmriprep/interfaces/__init__.py b/fmriprep/interfaces/__init__.py
index 32cdf5d1f..c8d29c87a 100644
--- a/fmriprep/interfaces/__init__.py
+++ b/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
diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py
index f02dae64c..50bf75786 100755
--- a/fmriprep/workflows/bold/base.py
+++ b/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):

From e5a6dee92eb9197313caaeb44b444b1bb42af51b Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Fri, 2 Nov 2018 12:17:01 -0400
Subject: [PATCH 09/34] Merge phase1phase2 into a single 4D file and send it to
 the DataSink

---
 fmriprep/workflows/fieldmap/phdiff.py | 28 +++++++++++++++------------
 1 file changed, 16 insertions(+), 12 deletions(-)

diff --git a/fmriprep/workflows/fieldmap/phdiff.py b/fmriprep/workflows/fieldmap/phdiff.py
index 158ed3e8b..c1ecbbb3b 100644
--- a/fmriprep/workflows/fieldmap/phdiff.py
+++ b/fmriprep/workflows/fieldmap/phdiff.py
@@ -75,17 +75,6 @@ def init_phdiff_wf(omp_nthreads, phasetype="phasediff", name='phdiff_wf'):
     def _pick1st(inlist):
         return inlist[0]
 
-    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
@@ -120,6 +109,22 @@ def _pick1st(inlist):
     #                     name='ComputeFieldmapFUGUE')
     # rsec2hz (divide by 2pi)
 
+    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')
+        workflow.connect(inputnode, 'phasediff', ds_fmap_mask, 'source_file')
+
+    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')
+        # Merge phase1 and phase2
+        phamrg = pe.Node(IntraModalMerge(), name='phamrg')
+        workflow.connect(inputnode, 'phasediff', phamrg, 'in_files')
+        workflow.connect(phamrg, 'out_file', ds_fmap_mask, 'source_file')
+
     workflow.connect([
         (inputnode, meta, [('phasediff', 'in_file')]),
         (inputnode, magmrg, [('magnitude', 'in_files')]),
@@ -138,7 +143,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')]),
     ])
 

From d89c776ec5ab42317341fd0fb089c46effb08a5d Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Fri, 2 Nov 2018 12:34:50 -0400
Subject: [PATCH 10/34] Send derived phasediff image to datasink

---
 fmriprep/interfaces/fmap.py           | 23 ++++++++++++++++++-----
 fmriprep/workflows/fieldmap/phdiff.py |  6 ++----
 2 files changed, 20 insertions(+), 9 deletions(-)

diff --git a/fmriprep/interfaces/fmap.py b/fmriprep/interfaces/fmap.py
index 143e8ef26..07b20273f 100644
--- a/fmriprep/interfaces/fmap.py
+++ b/fmriprep/interfaces/fmap.py
@@ -218,18 +218,25 @@ class Phases2FieldmapInputSpec(BaseInterfaceInputSpec):
     metadata = traits.List(mandatory=True, desc='BIDS metadata dictionaries')
 
 
+class Phases2FieldmapOutputSpec(TraitedSpec):
+    out_file = File(desc='the output fieldmap')
+    derived_phasediff = File(desc='the calculated phasediff (wrapped)')
+
+
 class Phases2Fieldmap(SimpleInterface):
     """
     Convert a phase difference map into a fieldmap in Hz
     """
     input_spec = Phases2FieldmapInputSpec
-    output_spec = Phasediff2FieldmapOutputSpec
+    output_spec = Phases2FieldmapOutputSpec
 
     def _run_interface(self, runtime):
-        self._results['out_file'] = phases2fmap(
+        phasediff_file, out_file = phases2fmap(
             self.inputs.in_file,
             _delta_te_from_two_phases(self.inputs.metadata),
             newpath=runtime.cwd)
+        self._results['out_file'] = out_file
+        self._results['derived_phasediff'] = phasediff_file
         return runtime
 
 
@@ -544,6 +551,7 @@ def phases2fmap(in_files, delta_te, newpath=None):
 
         \Delta B_0 (\text{Hz}) = \frac{\Delta \Theta}{2\pi \Delta\text{TE}}
 
+    Returns the calculated phasediff image and the fieldmap image
     """
     import math
     import numpy as np
@@ -552,14 +560,19 @@ def phases2fmap(in_files, delta_te, newpath=None):
     #  GYROMAG_RATIO_H_PROTON_MHZ = 42.576
 
     out_file = fname_presuffix(in_files[0], suffix='_fmap', newpath=newpath)
+    phasediff_file = fname_presuffix(in_files[0], suffix='_phasediff',
+                                     newpath=newpath)
     image0 = nb.load(in_files[0])
     image1 = nb.load(in_files[1])
-    data = image1.get_data() - image0.get_data()
-    data = (data / (2. * math.pi * delta_te))
+    phasediff = image1.get_data() - image0.get_data()
+    phasediff_nii = nb.Nifti1Image(phasediff, image0.affine, image0.header)
+    phasediff_nii.set_data_dtype(np.float32)
+    phasediff_nii.to_filename(phasediff_file)
+    data = (phasediff / (2. * math.pi * delta_te))
     nii = nb.Nifti1Image(data, image0.affine, image0.header)
     nii.set_data_dtype(np.float32)
     nii.to_filename(out_file)
-    return out_file
+    return phasediff_file, out_file
 
 
 def _delta_te(in_values, te1=None, te2=None):
diff --git a/fmriprep/workflows/fieldmap/phdiff.py b/fmriprep/workflows/fieldmap/phdiff.py
index c1ecbbb3b..e1ccf0ea7 100644
--- a/fmriprep/workflows/fieldmap/phdiff.py
+++ b/fmriprep/workflows/fieldmap/phdiff.py
@@ -120,10 +120,8 @@ def _pick1st(inlist):
         meta = pe.MapNode(ReadSidecarJSON(), name='meta', mem_gb=0.01,
                           run_without_submitting=True, iterfield=['in_file'])
         compfmap = pe.Node(Phases2Fieldmap(), name='compfmap')
-        # Merge phase1 and phase2
-        phamrg = pe.Node(IntraModalMerge(), name='phamrg')
-        workflow.connect(inputnode, 'phasediff', phamrg, 'in_files')
-        workflow.connect(phamrg, 'out_file', ds_fmap_mask, 'source_file')
+        workflow.connect(compfmap, 'derived_phasediff', ds_fmap_mask,
+                         'source_file')
 
     workflow.connect([
         (inputnode, meta, [('phasediff', 'in_file')]),

From f6d4aa71d926270f95a32552e58bf843288f0bfa Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Fri, 2 Nov 2018 12:53:57 -0400
Subject: [PATCH 11/34] Take phasediff from compfmap, not inputnode

---
 fmriprep/workflows/fieldmap/phdiff.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/fmriprep/workflows/fieldmap/phdiff.py b/fmriprep/workflows/fieldmap/phdiff.py
index e1ccf0ea7..108c8275c 100644
--- a/fmriprep/workflows/fieldmap/phdiff.py
+++ b/fmriprep/workflows/fieldmap/phdiff.py
@@ -115,6 +115,7 @@ def _pick1st(inlist):
                        run_without_submitting=True)
         compfmap = pe.Node(Phasediff2Fieldmap(), name='compfmap')
         workflow.connect(inputnode, 'phasediff', ds_fmap_mask, 'source_file')
+        workflow.connect(inputnode, 'phasediff', pha2rads, 'in_file')
 
     elif phasetype == "phase":
         meta = pe.MapNode(ReadSidecarJSON(), name='meta', mem_gb=0.01,
@@ -122,6 +123,7 @@ def _pick1st(inlist):
         compfmap = pe.Node(Phases2Fieldmap(), name='compfmap')
         workflow.connect(compfmap, 'derived_phasediff', ds_fmap_mask,
                          'source_file')
+        workflow.connect(compfmap, 'derived_phasediff', pha2rads, 'in_file')
 
     workflow.connect([
         (inputnode, meta, [('phasediff', 'in_file')]),
@@ -130,7 +132,7 @@ def _pick1st(inlist):
         (n4, prelude, [('output_image', 'magnitude_file')]),
         (n4, bet, [('output_image', 'in_file')]),
         (bet, prelude, [('mask_file', 'mask_file')]),
-        (inputnode, pha2rads, [('phasediff', 'in_file')]),
+
         (meta, compfmap, [('out_dict', 'metadata')]),
         (pha2rads, prelude, [('out', 'phase_file')]),
         (prelude, denoise, [('unwrapped_phase_file', 'in_file')]),

From 649dafe08759e0b90ad83fac68e991016161e567 Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Fri, 2 Nov 2018 14:03:47 -0400
Subject: [PATCH 12/34] Refactored interfaces

---
 fmriprep/interfaces/__init__.py       | 10 ++-
 fmriprep/interfaces/fmap.py           | 92 +++++++++++++--------------
 fmriprep/workflows/fieldmap/phdiff.py | 33 ++++++----
 3 files changed, 73 insertions(+), 62 deletions(-)

diff --git a/fmriprep/interfaces/__init__.py b/fmriprep/interfaces/__init__.py
index c8d29c87a..01ae821bf 100644
--- a/fmriprep/interfaces/__init__.py
+++ b/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,11 @@
 )
 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, Phases2Fieldmap
+from .utils import (TPM2ROI, AddTPMs, AddTSVHeader, ConcatAffines,
+                    JoinTSVColumns)
+from .fmap import (
+    FieldEnhance, FieldToRadS, FieldToHz, Phasediff2Fieldmap, Phases2Phasediff,
+    PhasesMetadata2PhasediffMetadata)
 from .confounds import GatherConfounds, ICAConfounds, FMRISummary
 from .itk import MCFLIRT2ITK, MultiApplyTransforms
 from .multiecho import FirstEcho
diff --git a/fmriprep/interfaces/fmap.py b/fmriprep/interfaces/fmap.py
index 07b20273f..98470e3ef 100644
--- a/fmriprep/interfaces/fmap.py
+++ b/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, InputMultiPath, traits)
 
 LOGGER = logging.getLogger('nipype.interface')
 
@@ -213,30 +213,46 @@ 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 Phases2PhasediffInputSpec(BaseInterfaceInputSpec):
+    in_files = InputMultiPath(File(exists=True), mandatory=True,
+                              desc='list of phase1, phase2 files')
 
 
-class Phases2FieldmapOutputSpec(TraitedSpec):
-    out_file = File(desc='the output fieldmap')
-    derived_phasediff = File(desc='the calculated phasediff (wrapped)')
+class Phases2PhasediffOutputSpec(TraitedSpec):
+    out_file = File(desc='the output phasediff')
 
 
-class Phases2Fieldmap(SimpleInterface):
+class Phases2Phasediff(SimpleInterface):
     """
-    Convert a phase difference map into a fieldmap in Hz
+    Convert a phase1, phase2 into a difference map
     """
-    input_spec = Phases2FieldmapInputSpec
-    output_spec = Phases2FieldmapOutputSpec
+    input_spec = Phases2PhasediffInputSpec
+    output_spec = Phases2PhasediffOutputSpec
 
     def _run_interface(self, runtime):
-        phasediff_file, out_file = phases2fmap(
-            self.inputs.in_file,
-            _delta_te_from_two_phases(self.inputs.metadata),
-            newpath=runtime.cwd)
-        self._results['out_file'] = out_file
-        self._results['derived_phasediff'] = phasediff_file
+        self._results['out_file'] = phases2phasediff(self.inputs.in_files)
+        return runtime
+
+
+class PhasesMetadata2PhasediffMetadataInputSpec(BaseInterfaceInputSpec):
+    in_dicts = traits.List(traits.Dict, mandatory=True,
+                           desc='list of phase1, phase2 metadata dicts')
+
+
+class PhasesMetadata2PhasediffMetadataOutputSpec(TraitedSpec):
+    out_dict = traits.Dict(desc='the phasediff metadata')
+
+
+class PhasesMetadata2PhasediffMetadata(SimpleInterface):
+    """
+    Convert phase1, phase2 metadata into phase difference metadata
+    """
+    input_spec = PhasesMetadata2PhasediffMetadataInputSpec
+    output_spec = PhasesMetadata2PhasediffMetadataOutputSpec
+
+    def _run_interface(self, runtime):
+        self._results['out_dict'] = phases_metadata_to_phasediff_metadata(
+                                        self.inputs.in_dicts)
         return runtime
 
 
@@ -534,32 +550,11 @@ def phdiff2fmap(in_file, delta_te, newpath=None):
     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}}
-
-    Returns the calculated phasediff image and the fieldmap image
-    """
-    import math
+def phases2phasediff(in_files, newpath=None):
+    """calculates a phasediff from two phase images"""
     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_files[0], suffix='_fmap', newpath=newpath)
     phasediff_file = fname_presuffix(in_files[0], suffix='_phasediff',
                                      newpath=newpath)
     image0 = nb.load(in_files[0])
@@ -568,11 +563,7 @@ def phases2fmap(in_files, delta_te, newpath=None):
     phasediff_nii = nb.Nifti1Image(phasediff, image0.affine, image0.header)
     phasediff_nii.set_data_dtype(np.float32)
     phasediff_nii.to_filename(phasediff_file)
-    data = (phasediff / (2. * math.pi * delta_te))
-    nii = nb.Nifti1Image(data, image0.affine, image0.header)
-    nii.set_data_dtype(np.float32)
-    nii.to_filename(out_file)
-    return phasediff_file, out_file
+    return phasediff_file
 
 
 def _delta_te(in_values, te1=None, te2=None):
@@ -612,14 +603,16 @@ def _delta_te(in_values, te1=None, te2=None):
     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"""
+def phases_metadata_to_phasediff_metadata(in_dicts, te1=None, te2=None):
+    """Combines two BIDS metadata dicts into a phasediff-like metadata dict"""
 
+    reference_dict = {"EchoTime": None}
     if isinstance(in_dicts, list):
         te2_value, te1_value = in_dicts
         if isinstance(te1_value, list):
             te1_dict = te1[1]
             if isinstance(te1_dict, dict):
+                reference_dict = te1_dict
                 te1 = te1_dict.get('EchoTime')
         if isinstance(te2_value, list):
             te2_dict = te2[1]
@@ -639,4 +632,7 @@ def _delta_te_from_two_phases(in_dicts, te1=None, te2=None):
             'EchoTime metadata field not found for phase2. Please consult the '
             'BIDS specification.')
 
-    return float(te2) - float(te1)
+    del reference_dict['EchoTime']
+    reference_dict["EchoTime1"] = te1
+    reference_dict["EchoTime2"] = te2
+    return reference_dict
diff --git a/fmriprep/workflows/fieldmap/phdiff.py b/fmriprep/workflows/fieldmap/phdiff.py
index 108c8275c..891467707 100644
--- a/fmriprep/workflows/fieldmap/phdiff.py
+++ b/fmriprep/workflows/fieldmap/phdiff.py
@@ -27,7 +27,7 @@
 from ...engine import Workflow
 from ...interfaces import (
     ReadSidecarJSON, IntraModalMerge, DerivativesDataSink,
-    Phasediff2Fieldmap, Phases2Fieldmap
+    Phasediff2Fieldmap, Phases2Phasediff, PhasesMetadata2PhasediffMetadata
 )
 
 
@@ -103,6 +103,8 @@ 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)
     # pre_fugue = pe.Node(fsl.FUGUE(save_fmap=True),
@@ -113,27 +115,36 @@ def _pick1st(inlist):
         # Read phasediff echo times
         meta = pe.Node(ReadSidecarJSON(), name='meta', mem_gb=0.01,
                        run_without_submitting=True)
-        compfmap = pe.Node(Phasediff2Fieldmap(), name='compfmap')
-        workflow.connect(inputnode, 'phasediff', ds_fmap_mask, 'source_file')
-        workflow.connect(inputnode, 'phasediff', pha2rads, 'in_file')
+        workflow.connect([
+
+            (meta, compfmap, [('out_dict', 'metadata')]),
+            (inputnode, pha2rads, [('phasediff', 'in_file')]),
+            (inputnode, ds_fmap_mask, [('phasediff', 'source_file')]),
+        ])
 
     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')
-        workflow.connect(compfmap, 'derived_phasediff', ds_fmap_mask,
-                         'source_file')
-        workflow.connect(compfmap, 'derived_phasediff', pha2rads, 'in_file')
+        phasediff_meta = pe.Node(PhasesMetadata2PhasediffMetadata(),
+                                 name='meta', mem_gb=0.01,
+                                 run_without_submitting=True)
+        diff_phases = pe.Node(Phases2Phasediff(), name='diff_phases')
+
+        workflow.connect([
+            (inputnode, meta, [('phasediff', 'in_file')]),
+            (meta, phasediff_meta, [('out_dice', 'in_file')]),
+            (phasediff_meta, compfmap, [('out_dict', 'metadata')]),
+            (inputnode, diff_phases, [('phasediff', 'in_files')]),
+            (diff_phases, pha2rads, [('out_file', 'in_file')]),
+            (diff_phases, 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')]),
-
-        (meta, compfmap, [('out_dict', 'metadata')]),
         (pha2rads, prelude, [('out', 'phase_file')]),
         (prelude, denoise, [('unwrapped_phase_file', 'in_file')]),
         (denoise, demean, [('out_file', 'in_file')]),

From 8558aa52e49ac220b05b0fcb2c166ef7a999da1e Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Fri, 2 Nov 2018 15:55:30 -0400
Subject: [PATCH 13/34] fix typo and node name

---
 fmriprep/interfaces/fmap.py           | 23 +++++++++--------------
 fmriprep/workflows/fieldmap/phdiff.py |  4 ++--
 2 files changed, 11 insertions(+), 16 deletions(-)

diff --git a/fmriprep/interfaces/fmap.py b/fmriprep/interfaces/fmap.py
index 98470e3ef..096406a4a 100644
--- a/fmriprep/interfaces/fmap.py
+++ b/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, InputMultiPath, traits)
+    SimpleInterface, InputMultiPath)
 
 LOGGER = logging.getLogger('nipype.interface')
 
@@ -559,8 +559,9 @@ def phases2phasediff(in_files, newpath=None):
                                      newpath=newpath)
     image0 = nb.load(in_files[0])
     image1 = nb.load(in_files[1])
-    phasediff = image1.get_data() - image0.get_data()
-    phasediff_nii = nb.Nifti1Image(phasediff, image0.affine, image0.header)
+    phasediff = image1.get_data().astype(np.float) \
+                - image0.get_data().astype(np.float)
+    phasediff_nii = nb.Nifti1Image(phasediff, image0.affine)
     phasediff_nii.set_data_dtype(np.float32)
     phasediff_nii.to_filename(phasediff_file)
     return phasediff_file
@@ -605,19 +606,13 @@ def _delta_te(in_values, te1=None, te2=None):
 
 def phases_metadata_to_phasediff_metadata(in_dicts, te1=None, te2=None):
     """Combines two BIDS metadata dicts into a phasediff-like metadata dict"""
+    from copy import deepcopy
 
     reference_dict = {"EchoTime": None}
-    if isinstance(in_dicts, list):
-        te2_value, te1_value = in_dicts
-        if isinstance(te1_value, list):
-            te1_dict = te1[1]
-            if isinstance(te1_dict, dict):
-                reference_dict = te1_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')
+    phase1_dict, phase2_dict = in_dicts
+    reference_dict = deepcopy(phase1_dict)
+    te1 = phase1_dict.get('EchoTime')
+    te2 = phase2_dict.get('EchoTime')
 
     # For convienience if both are missing we should give one error about them
     if te1 is None and te2 is None:
diff --git a/fmriprep/workflows/fieldmap/phdiff.py b/fmriprep/workflows/fieldmap/phdiff.py
index 891467707..eb21f03db 100644
--- a/fmriprep/workflows/fieldmap/phdiff.py
+++ b/fmriprep/workflows/fieldmap/phdiff.py
@@ -126,13 +126,13 @@ def _pick1st(inlist):
         meta = pe.MapNode(ReadSidecarJSON(), name='meta', mem_gb=0.01,
                           run_without_submitting=True, iterfield=['in_file'])
         phasediff_meta = pe.Node(PhasesMetadata2PhasediffMetadata(),
-                                 name='meta', mem_gb=0.01,
+                                 name='phasediff_meta', mem_gb=0.01,
                                  run_without_submitting=True)
         diff_phases = pe.Node(Phases2Phasediff(), name='diff_phases')
 
         workflow.connect([
             (inputnode, meta, [('phasediff', 'in_file')]),
-            (meta, phasediff_meta, [('out_dice', 'in_file')]),
+            (meta, phasediff_meta, [('out_dict', 'in_dicts')]),
             (phasediff_meta, compfmap, [('out_dict', 'metadata')]),
             (inputnode, diff_phases, [('phasediff', 'in_files')]),
             (diff_phases, pha2rads, [('out_file', 'in_file')]),

From 1d12256262df842c34d10cd991d127c08dea8367 Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Sat, 3 Nov 2018 12:30:42 -0400
Subject: [PATCH 14/34] Calculate input to prelude correctly

---
 fmriprep/interfaces/__init__.py       |   3 +-
 fmriprep/interfaces/fmap.py           | 129 ++++++++++++--------------
 fmriprep/workflows/fieldmap/phdiff.py |  25 +++--
 3 files changed, 74 insertions(+), 83 deletions(-)

diff --git a/fmriprep/interfaces/__init__.py b/fmriprep/interfaces/__init__.py
index 01ae821bf..7c332b2a2 100644
--- a/fmriprep/interfaces/__init__.py
+++ b/fmriprep/interfaces/__init__.py
@@ -17,8 +17,7 @@
 from .utils import (TPM2ROI, AddTPMs, AddTSVHeader, ConcatAffines,
                     JoinTSVColumns)
 from .fmap import (
-    FieldEnhance, FieldToRadS, FieldToHz, Phasediff2Fieldmap, Phases2Phasediff,
-    PhasesMetadata2PhasediffMetadata)
+    FieldEnhance, FieldToRadS, FieldToHz, Phasediff2Fieldmap, Phases2Fieldmap)
 from .confounds import GatherConfounds, ICAConfounds, FMRISummary
 from .itk import MCFLIRT2ITK, MultiApplyTransforms
 from .multiecho import FirstEcho
diff --git a/fmriprep/interfaces/fmap.py b/fmriprep/interfaces/fmap.py
index 096406a4a..9eb2eef6f 100644
--- a/fmriprep/interfaces/fmap.py
+++ b/fmriprep/interfaces/fmap.py
@@ -213,46 +213,35 @@ def _run_interface(self, runtime):
         return runtime
 
 
-class Phases2PhasediffInputSpec(BaseInterfaceInputSpec):
-    in_files = InputMultiPath(File(exists=True), mandatory=True,
-                              desc='list of phase1, phase2 files')
+class Phases2FieldmapInputSpec(BaseInterfaceInputSpec):
+    phase_files = InputMultiPath(File(exists=True), mandatory=True,
+                                 desc='list of phase1, phase2 files')
+    magnitude_files = InputMultiPath(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 Phases2PhasediffOutputSpec(TraitedSpec):
-    out_file = File(desc='the output phasediff')
+class Phases2FieldmapOutputSpec(TraitedSpec):
+    out_file = File(desc='the output fieldmap')
+    phasediff_metadata = traits.Dict(desc='the phasediff metadata')
 
 
-class Phases2Phasediff(SimpleInterface):
+class Phases2Fieldmap(SimpleInterface):
     """
     Convert a phase1, phase2 into a difference map
     """
-    input_spec = Phases2PhasediffInputSpec
-    output_spec = Phases2PhasediffOutputSpec
-
-    def _run_interface(self, runtime):
-        self._results['out_file'] = phases2phasediff(self.inputs.in_files)
-        return runtime
-
-
-class PhasesMetadata2PhasediffMetadataInputSpec(BaseInterfaceInputSpec):
-    in_dicts = traits.List(traits.Dict, mandatory=True,
-                           desc='list of phase1, phase2 metadata dicts')
-
-
-class PhasesMetadata2PhasediffMetadataOutputSpec(TraitedSpec):
-    out_dict = traits.Dict(desc='the phasediff metadata')
-
-
-class PhasesMetadata2PhasediffMetadata(SimpleInterface):
-    """
-    Convert phase1, phase2 metadata into phase difference metadata
-    """
-    input_spec = PhasesMetadata2PhasediffMetadataInputSpec
-    output_spec = PhasesMetadata2PhasediffMetadataOutputSpec
+    input_spec = Phases2FieldmapInputSpec
+    output_spec = Phases2FieldmapOutputSpec
 
     def _run_interface(self, runtime):
-        self._results['out_dict'] = phases_metadata_to_phasediff_metadata(
-                                        self.inputs.in_dicts)
+        # Get the echo times
+        fmap_file, merged_metadata = phases2fmap(
+            self.inputs.phase_files, self.inputs.magnitude_files,
+            self.inputs.metadatas
+        )
+        self._results['phasediff_metadata'] = merged_metadata
+        self._results['out_file'] = fmap_file
         return runtime
 
 
@@ -550,20 +539,53 @@ def phdiff2fmap(in_file, delta_te, newpath=None):
     return out_file
 
 
-def phases2phasediff(in_files, newpath=None):
-    """calculates a phasediff from two phase images"""
+def phases2fmap(phase_files, magnitude_files, metadatas, newpath=None):
+    """calculates a phasediff from two phase images. Assumes monopolar
+    readout and that the second image in the list corresponds to the
+    longer echo time"""
     import numpy as np
     import nibabel as nb
     from nipype.utils.filemanip import fname_presuffix
-    phasediff_file = fname_presuffix(in_files[0], suffix='_phasediff',
+    from copy import deepcopy
+
+    phasediff_file = fname_presuffix(phase_files[0], suffix='_phasediff',
                                      newpath=newpath)
-    image0 = nb.load(in_files[0])
-    image1 = nb.load(in_files[1])
-    phasediff = image1.get_data().astype(np.float) \
-                - image0.get_data().astype(np.float)
-    phasediff_nii = nb.Nifti1Image(phasediff, image0.affine)
+    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))
+
+    magnitude_image = magnitude_files[short_echo_index]
+    short_phase_image = phase_files[short_echo_index]
+    long_phase_image = phase_files[long_echo_index]
+
+    mag_img = nb.load(magnitude_image)
+    mag = mag_img.get_data().astype(np.float)
+    image0 = nb.load(short_phase_image)
+    phase0 = image0.get_data().astype(np.float)
+    image1 = nb.load(long_phase_image)
+    phase1 = image1.get_data().astype(np.float)
+
+    # Calculate fieldmaps
+    rad0 = np.pi * phase0 / 2048 - np.pi
+    rad1 = np.pi * phase1 / 2048 - np.pi
+    a = mag * np.cos(rad0)
+    b = mag * np.sin(rad0)
+    c = mag * np.cos(rad1)
+    d = mag * np.sin(rad1)
+    fmap = -np.atan2(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
 
 
@@ -602,32 +624,3 @@ def _delta_te(in_values, te1=None, te2=None):
             'specification.')
 
     return abs(float(te2) - float(te1))
-
-
-def phases_metadata_to_phasediff_metadata(in_dicts, te1=None, te2=None):
-    """Combines two BIDS metadata dicts into a phasediff-like metadata dict"""
-    from copy import deepcopy
-
-    reference_dict = {"EchoTime": None}
-    phase1_dict, phase2_dict = in_dicts
-    reference_dict = deepcopy(phase1_dict)
-    te1 = phase1_dict.get('EchoTime')
-    te2 = phase2_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.')
-
-    del reference_dict['EchoTime']
-    reference_dict["EchoTime1"] = te1
-    reference_dict["EchoTime2"] = te2
-    return reference_dict
diff --git a/fmriprep/workflows/fieldmap/phdiff.py b/fmriprep/workflows/fieldmap/phdiff.py
index eb21f03db..698d79ee9 100644
--- a/fmriprep/workflows/fieldmap/phdiff.py
+++ b/fmriprep/workflows/fieldmap/phdiff.py
@@ -27,7 +27,7 @@
 from ...engine import Workflow
 from ...interfaces import (
     ReadSidecarJSON, IntraModalMerge, DerivativesDataSink,
-    Phasediff2Fieldmap, Phases2Phasediff, PhasesMetadata2PhasediffMetadata
+    Phasediff2Fieldmap, Phases2Fieldmap
 )
 
 
@@ -75,8 +75,6 @@ def init_phdiff_wf(omp_nthreads, phasetype="phasediff", name='phdiff_wf'):
     def _pick1st(inlist):
         return inlist[0]
 
-    pha2rads = pe.Node(niu.Function(function=siemens2rads), name='pha2rads')
-
     # Merge input magnitude images
     magmrg = pe.Node(IntraModalMerge(), name='magmrg')
 
@@ -112,6 +110,9 @@ def _pick1st(inlist):
     # 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)
@@ -119,24 +120,23 @@ def _pick1st(inlist):
 
             (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":
         meta = pe.MapNode(ReadSidecarJSON(), name='meta', mem_gb=0.01,
                           run_without_submitting=True, iterfield=['in_file'])
-        phasediff_meta = pe.Node(PhasesMetadata2PhasediffMetadata(),
-                                 name='phasediff_meta', mem_gb=0.01,
-                                 run_without_submitting=True)
-        diff_phases = pe.Node(Phases2Phasediff(), name='diff_phases')
+        phases2fmap = pe.Node(Phases2Fieldmap(), name='phases2fmap')
 
         workflow.connect([
             (inputnode, meta, [('phasediff', 'in_file')]),
-            (meta, phasediff_meta, [('out_dict', 'in_dicts')]),
-            (phasediff_meta, compfmap, [('out_dict', 'metadata')]),
-            (inputnode, diff_phases, [('phasediff', 'in_files')]),
-            (diff_phases, pha2rads, [('out_file', 'in_file')]),
-            (diff_phases, ds_fmap_mask, [('out_file', 'source_file')])
+            (meta, phases2fmap, [('out_dict', 'metadatas')]),
+            (inputnode, phases2fmap, [('phasediff', 'phase_files')]),
+            (inputnode, phases2fmap, [('magnitude', 'magnitude_files')]),
+            (phases2fmap, prelude), [('out_file', 'phase_file')],
+            (phases2fmap, compfmap, [('phasediff_metadata', 'metadata')]),
+            (phases2fmap, ds_fmap_mask, [('out_file', 'source_file')])
         ])
 
     workflow.connect([
@@ -145,7 +145,6 @@ def _pick1st(inlist):
         (n4, prelude, [('output_image', 'magnitude_file')]),
         (n4, bet, [('output_image', 'in_file')]),
         (bet, prelude, [('mask_file', 'mask_file')]),
-        (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')]),

From 70b34131dfdf3789d989b1ecd6e0d8bcc06c56d3 Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Sun, 4 Nov 2018 09:02:12 -0500
Subject: [PATCH 15/34] phase1, phase2 worksgit status

---
 fmriprep/interfaces/fmap.py           |  6 +++---
 fmriprep/workflows/fieldmap/phdiff.py | 10 +++-------
 2 files changed, 6 insertions(+), 10 deletions(-)

diff --git a/fmriprep/interfaces/fmap.py b/fmriprep/interfaces/fmap.py
index 9eb2eef6f..fa4025473 100644
--- a/fmriprep/interfaces/fmap.py
+++ b/fmriprep/interfaces/fmap.py
@@ -575,18 +575,18 @@ def phases2fmap(phase_files, magnitude_files, metadatas, newpath=None):
     b = mag * np.sin(rad0)
     c = mag * np.cos(rad1)
     d = mag * np.sin(rad1)
-    fmap = -np.atan2(b*c-a*d, a*c+b*d)
+    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]]
+    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
+    return phasediff_file, merged_metadata
 
 
 def _delta_te(in_values, te1=None, te2=None):
diff --git a/fmriprep/workflows/fieldmap/phdiff.py b/fmriprep/workflows/fieldmap/phdiff.py
index 698d79ee9..b6a545e7a 100644
--- a/fmriprep/workflows/fieldmap/phdiff.py
+++ b/fmriprep/workflows/fieldmap/phdiff.py
@@ -72,9 +72,6 @@ def init_phdiff_wf(omp_nthreads, phasetype="phasediff", name='phdiff_wf'):
     outputnode = pe.Node(niu.IdentityInterface(
         fields=['fmap', 'fmap_ref', 'fmap_mask']), name='outputnode')
 
-    def _pick1st(inlist):
-        return inlist[0]
-
     # Merge input magnitude images
     magmrg = pe.Node(IntraModalMerge(), name='magmrg')
 
@@ -117,7 +114,6 @@ def _pick1st(inlist):
         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')]),
@@ -125,21 +121,21 @@ def _pick1st(inlist):
         ])
 
     elif phasetype == "phase":
+        # 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([
-            (inputnode, meta, [('phasediff', 'in_file')]),
             (meta, phases2fmap, [('out_dict', 'metadatas')]),
             (inputnode, phases2fmap, [('phasediff', 'phase_files')]),
             (inputnode, phases2fmap, [('magnitude', 'magnitude_files')]),
-            (phases2fmap, prelude), [('out_file', 'phase_file')],
+            (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')]),

From 191e1559096ea67b4b39ef7cae3b2520d54a47a2 Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Sun, 4 Nov 2018 09:05:48 -0500
Subject: [PATCH 16/34] Update zendodo json

---
 .zenodo.json | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/.zenodo.json b/.zenodo.json
index b1c31a81b..4dbad4bf4 100644
--- a/.zenodo.json
+++ b/.zenodo.json
@@ -97,6 +97,11 @@
       "name": "Gorgolewski, Krzysztof J.",
       "affiliation": "Department of Psychology, Stanford University",
       "orcid": "0000-0003-3321-7583"
+    },
+    {
+      "name": "Cieslak, Matthew",
+      "affiliation": "Department of Neuropsychiatry, University of Pennsylvania",
+      "orcid": "0000-0002-1931-4734"
     }
   ],
   "keywords": [

From 880d04d48d0b0c87002fdbd75b60ba9b8853bc1b Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Thu, 8 Nov 2018 12:34:03 -0500
Subject: [PATCH 17/34] Changes suggested during PR review

---
 .zenodo.json                          |  10 +--
 fmriprep/interfaces/fmap.py           | 116 ++++++++++----------------
 fmriprep/workflows/fieldmap/phdiff.py |   2 +-
 3 files changed, 51 insertions(+), 77 deletions(-)

diff --git a/.zenodo.json b/.zenodo.json
index 4dbad4bf4..989d2ad1f 100644
--- a/.zenodo.json
+++ b/.zenodo.json
@@ -83,6 +83,11 @@
       "affiliation": "Department of Psychology, University of Oslo",
       "orcid": "0000-0002-8508-9088"
     },
+    {
+      "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",
@@ -97,11 +102,6 @@
       "name": "Gorgolewski, Krzysztof J.",
       "affiliation": "Department of Psychology, Stanford University",
       "orcid": "0000-0003-3321-7583"
-    },
-    {
-      "name": "Cieslak, Matthew",
-      "affiliation": "Department of Neuropsychiatry, University of Pennsylvania",
-      "orcid": "0000-0002-1931-4734"
     }
   ],
   "keywords": [
diff --git a/fmriprep/interfaces/fmap.py b/fmriprep/interfaces/fmap.py
index fa4025473..9aa950748 100644
--- a/fmriprep/interfaces/fmap.py
+++ b/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, InputMultiPath)
+from nipype.interfaces.base import (BaseInterfaceInputSpec, TraitedSpec, File, isdefined, traits,
+                                    SimpleInterface, InputMultiObject)
 
 LOGGER = logging.getLogger('nipype.interface')
 
@@ -32,13 +31,10 @@ class FieldEnhanceInputSpec(BaseInterfaceInputSpec):
     in_magnitude = File(exists=True, desc='input magnitude')
     unwrap = traits.Bool(False, usedefault=True, desc='run phase unwrap')
     despike = traits.Bool(True, usedefault=True, desc='run despike filter')
-    bspline_smooth = traits.Bool(True, usedefault=True,
-                                 desc='run 3D bspline smoother')
+    bspline_smooth = traits.Bool(True, usedefault=True, desc='run 3D bspline smoother')
     mask_erode = traits.Int(1, usedefault=True, desc='mask erosion iterations')
-    despike_threshold = traits.Float(0.2, usedefault=True,
-                                     desc='mask erosion iterations')
-    num_threads = traits.Int(1, usedefault=True, nohash=True,
-                             desc='number of jobs')
+    despike_threshold = traits.Float(0.2, usedefault=True, desc='mask erosion iterations')
+    num_threads = traits.Int(1, usedefault=True, nohash=True, desc='number of jobs')
 
 
 class FieldEnhanceOutputSpec(TraitedSpec):
@@ -71,12 +67,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)
+                struc = sim.iterate_structure(sim.generate_binary_structure(3, 2), 1)
                 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)  # pylint: disable=no-member
 
         self._results['out_file'] = fname_presuffix(
             self.inputs.in_file, suffix='_enh', newpath=runtime.cwd)
@@ -86,8 +79,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'])
@@ -97,8 +90,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()
 
@@ -112,8 +104,7 @@ def _run_interface(self, runtime):
 
             nslices = 0
             try:
-                errorslice = np.squeeze(
-                                np.argwhere(errormask.sum(0).sum(0) > 0))
+                errorslice = np.squeeze(np.argwhere(errormask.sum(0).sum(0) > 0))
                 nslices = errorslice[-1] - errorslice[0]
             except IndexError:  # mask is empty, do not refine
                 pass
@@ -121,12 +112,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()
 
@@ -135,8 +125,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
 
@@ -207,19 +197,15 @@ 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 = InputMultiPath(File(exists=True), mandatory=True,
-                                 desc='list of phase1, phase2 files')
-    magnitude_files = InputMultiPath(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')
+    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):
@@ -236,10 +222,7 @@ class Phases2Fieldmap(SimpleInterface):
 
     def _run_interface(self, runtime):
         # Get the echo times
-        fmap_file, merged_metadata = phases2fmap(
-            self.inputs.phase_files, self.inputs.magnitude_files,
-            self.inputs.metadatas
-        )
+        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
@@ -271,8 +254,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
 
@@ -290,16 +272,15 @@ def _unwrap(fmap_data, mag_file, mask=None):
 
     nb.Nifti1Image(fmap_data, magnii.affine).to_filename('fmap_rad.nii.gz')
     nb.Nifti1Image(mask, magnii.affine).to_filename('fmap_mask.nii.gz')
-    nb.Nifti1Image(magnii.get_data(), magnii.affine
-                   ).to_filename('fmap_mag.nii.gz')
+    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)
+    unwrapped = nb.load(res.outputs.unwrapped_phase_file).get_data() * (fmapmax / pi)
     return unwrapped
 
 
@@ -539,17 +520,15 @@ def phdiff2fmap(in_file, delta_te, newpath=None):
     return out_file
 
 
-def phases2fmap(phase_files, magnitude_files, metadatas, newpath=None):
-    """calculates a phasediff from two phase images. Assumes monopolar
-    readout and that the second image in the list corresponds to the
-    longer echo time"""
+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)
+    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()
@@ -557,25 +536,22 @@ def phases2fmap(phase_files, magnitude_files, metadatas, newpath=None):
     short_echo_index = echo_times.index(min(echo_times))
     long_echo_index = echo_times.index(max(echo_times))
 
-    magnitude_image = magnitude_files[short_echo_index]
     short_phase_image = phase_files[short_echo_index]
     long_phase_image = phase_files[long_echo_index]
 
-    mag_img = nb.load(magnitude_image)
-    mag = mag_img.get_data().astype(np.float)
     image0 = nb.load(short_phase_image)
-    phase0 = image0.get_data().astype(np.float)
+    phase0 = image0.get_fdata()
     image1 = nb.load(long_phase_image)
-    phase1 = image1.get_data().astype(np.float)
+    phase1 = image1.get_fdata()
 
     # Calculate fieldmaps
     rad0 = np.pi * phase0 / 2048 - np.pi
     rad1 = np.pi * phase1 / 2048 - np.pi
-    a = mag * np.cos(rad0)
-    b = mag * np.sin(rad0)
-    c = mag * np.cos(rad1)
-    d = mag * np.sin(rad1)
-    fmap = -np.arctan2(b*c-a*d, a*c+b*d)
+    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)
@@ -615,12 +591,10 @@ def _delta_te(in_values, te1=None, te2=None):
         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))
diff --git a/fmriprep/workflows/fieldmap/phdiff.py b/fmriprep/workflows/fieldmap/phdiff.py
index b6a545e7a..b03f38d2d 100644
--- a/fmriprep/workflows/fieldmap/phdiff.py
+++ b/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
 
 """
 
@@ -128,7 +129,6 @@ def init_phdiff_wf(omp_nthreads, phasetype="phasediff", name='phdiff_wf'):
         workflow.connect([
             (meta, phases2fmap, [('out_dict', 'metadatas')]),
             (inputnode, phases2fmap, [('phasediff', 'phase_files')]),
-            (inputnode, phases2fmap, [('magnitude', 'magnitude_files')]),
             (phases2fmap, prelude, [('out_file', 'phase_file')]),
             (phases2fmap, compfmap, [('phasediff_metadata', 'metadata')]),
             (phases2fmap, ds_fmap_mask, [('out_file', 'source_file')])

From 9f9284ec9808bcc7aa37bececeb8a64f9fbc8077 Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Thu, 8 Nov 2018 12:40:23 -0500
Subject: [PATCH 18/34] line was too long

---
 fmriprep/interfaces/fmap.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/fmriprep/interfaces/fmap.py b/fmriprep/interfaces/fmap.py
index 9aa950748..563d674f3 100644
--- a/fmriprep/interfaces/fmap.py
+++ b/fmriprep/interfaces/fmap.py
@@ -68,8 +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)

From 0712c5e0e6b237e8b9b2b9d3aaaad415b0cff66d Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Thu, 8 Nov 2018 13:08:49 -0500
Subject: [PATCH 19/34] Add missing bracket

---
 .zenodo.json | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.zenodo.json b/.zenodo.json
index a456927d1..e3afdb33a 100644
--- a/.zenodo.json
+++ b/.zenodo.json
@@ -88,6 +88,7 @@
       "affiliation": "Department of Neuropsychiatry, University of Pennsylvania",
       "orcid": "0000-0002-1931-4734"
     },
+    {
       "name": "Liem, Franz",
       "affiliation": "University of Zurich",
       "orcid": "0000-0003-0646-4810"

From 09316a927f29a65b8612a659d04c31ebeddddf95 Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Thu, 8 Nov 2018 14:16:30 -0500
Subject: [PATCH 20/34] Make scaling more general

---
 fmriprep/interfaces/fmap.py | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/fmriprep/interfaces/fmap.py b/fmriprep/interfaces/fmap.py
index 563d674f3..7cad920b0 100644
--- a/fmriprep/interfaces/fmap.py
+++ b/fmriprep/interfaces/fmap.py
@@ -545,9 +545,15 @@ def phases2fmap(phase_files, metadatas, newpath=None):
     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 = np.pi * phase0 / 2048 - np.pi
-    rad1 = np.pi * phase1 / 2048 - np.pi
+    rad0 = rescale_image(phase0)
+    rad1 = rescale_image(phase1)
     a = np.cos(rad0)
     b = np.sin(rad0)
     c = np.cos(rad1)

From f9ca876993c7788fdfb05df7960c1f405b4e0ff0 Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Thu, 8 Nov 2018 16:06:54 -0500
Subject: [PATCH 21/34] Credit original script author

---
 fmriprep/workflows/fieldmap/phdiff.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/fmriprep/workflows/fieldmap/phdiff.py b/fmriprep/workflows/fieldmap/phdiff.py
index b03f38d2d..0db86d809 100644
--- a/fmriprep/workflows/fieldmap/phdiff.py
+++ b/fmriprep/workflows/fieldmap/phdiff.py
@@ -122,6 +122,10 @@ def init_phdiff_wf(omp_nthreads, phasetype="phasediff", name='phdiff_wf'):
         ])
 
     elif phasetype == "phase":
+        workflow.__desc__ += """\
+The phase difference used for unwarping was calculated using a method developed
+by Mark Elliot (University of Pennsylvania).
+    """
         # Special case for phase1, phase2 images
         meta = pe.MapNode(ReadSidecarJSON(), name='meta', mem_gb=0.01,
                           run_without_submitting=True, iterfield=['in_file'])

From ada0d9dc28cc58a4027c70c872d3c70e754c9643 Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Tue, 20 Nov 2018 16:15:43 -0500
Subject: [PATCH 22/34] Undo formatting changes and add citation

---
 fmriprep/data/boilerplate.bib         | 10 +++++++
 fmriprep/interfaces/__init__.py       |  9 ++----
 fmriprep/interfaces/fmap.py           | 41 +++++++++++++--------------
 fmriprep/workflows/fieldmap/base.py   | 38 +++++++++----------------
 fmriprep/workflows/fieldmap/phdiff.py | 14 ++++-----
 5 files changed, 51 insertions(+), 61 deletions(-)

diff --git a/fmriprep/data/boilerplate.bib b/fmriprep/data/boilerplate.bib
index a14f5cd1d..efd2f8bd7 100644
--- a/fmriprep/data/boilerplate.bib
+++ b/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},
diff --git a/fmriprep/interfaces/__init__.py b/fmriprep/interfaces/__init__.py
index 7c332b2a2..c8d29c87a 100644
--- a/fmriprep/interfaces/__init__.py
+++ b/fmriprep/interfaces/__init__.py
@@ -2,8 +2,7 @@
 # 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
@@ -14,10 +13,8 @@
 )
 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, Phases2Fieldmap)
+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
diff --git a/fmriprep/interfaces/fmap.py b/fmriprep/interfaces/fmap.py
index 7cad920b0..ead71a875 100644
--- a/fmriprep/interfaces/fmap.py
+++ b/fmriprep/interfaces/fmap.py
@@ -19,8 +19,9 @@
 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, InputMultiObject)
+from nipype.interfaces.base import (
+    BaseInterfaceInputSpec, TraitedSpec, File, isdefined, traits,
+    SimpleInterface, InputMultiObject)
 
 LOGGER = logging.getLogger('nipype.interface')
 
@@ -68,9 +69,10 @@ 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)
+                    mask, struc,
+                    iterations=self.inputs.mask_erode
+                ).astype(np.uint8) # pylint: disable=no-member
 
         self._results['out_file'] = fname_presuffix(
             self.inputs.in_file, suffix='_enh', newpath=runtime.cwd)
@@ -80,8 +82,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'])
@@ -276,10 +278,9 @@ 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
@@ -312,9 +313,8 @@ 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.).
 
@@ -494,8 +494,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
@@ -595,13 +594,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))
diff --git a/fmriprep/workflows/fieldmap/base.py b/fmriprep/workflows/fieldmap/base.py
index e3d780925..e6f79d666 100644
--- a/fmriprep/workflows/fieldmap/base.py
+++ b/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**)
 
@@ -81,17 +81,13 @@ 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',
             },
         )
@@ -180,12 +176,10 @@ 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,
@@ -219,16 +213,11 @@ def init_sdc_wf(fmaps, bold_meta, omp_nthreads=1,
                                                phasetype=fmap['type'])
             # set inputs
             if fmap['type'] == 'phasediff':
-                fmap_estimator_wf.inputs.inputnode.phasediff = \
-                    fmap['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.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,
@@ -257,8 +246,7 @@ 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')]),
diff --git a/fmriprep/workflows/fieldmap/phdiff.py b/fmriprep/workflows/fieldmap/phdiff.py
index 0db86d809..6b5cc43ca 100644
--- a/fmriprep/workflows/fieldmap/phdiff.py
+++ b/fmriprep/workflows/fieldmap/phdiff.py
@@ -66,8 +66,7 @@ def init_phdiff_wf(omp_nthreads, phasetype="phasediff", name='phdiff_wf'):
 further improvements of HCP Pipelines [@hcppipelines].
 """
 
-    inputnode = pe.Node(niu.IdentityInterface(
-                        fields=['magnitude', 'phasediff']),
+    inputnode = pe.Node(niu.IdentityInterface(fields=['magnitude', 'phasediff']),
                         name='inputnode')
 
     outputnode = pe.Node(niu.IdentityInterface(
@@ -81,8 +80,7 @@ def init_phdiff_wf(omp_nthreads, phasetype="phasediff", name='phdiff_wf'):
                  name='n4', n_procs=omp_nthreads)
     bet = pe.Node(BETRPT(generate_report=True, frac=0.6, mask=True),
                   name='bet')
-    ds_fmap_mask = pe.Node(DerivativesDataSink(suffix='fmap_mask'),
-                           name='ds_report_fmap_mask',
+    ds_fmap_mask = pe.Node(DerivativesDataSink(suffix='fmap_mask'), name='ds_report_fmap_mask',
                            mem_gb=0.01, run_without_submitting=True)
     # uses mask from bet; outputs a mask
     # dilate = pe.Node(fsl.maths.MathsCommand(
@@ -91,8 +89,7 @@ def init_phdiff_wf(omp_nthreads, phasetype="phasediff", name='phdiff_wf'):
     # FSL PRELUDE will perform phase-unwrapping
     prelude = pe.Node(fsl.PRELUDE(), name='prelude')
 
-    denoise = pe.Node(fsl.SpatialFilter(operation='median',
-                                        kernel_shape='sphere',
+    denoise = pe.Node(fsl.SpatialFilter(operation='median', kernel_shape='sphere',
                                         kernel_size=3), name='denoise')
 
     demean = pe.Node(niu.Function(function=demean_image), name='demean')
@@ -103,8 +100,7 @@ def init_phdiff_wf(omp_nthreads, phasetype="phasediff", name='phdiff_wf'):
 
     # The phdiff2fmap interface is equivalent to:
     # rad2rsec (using rads2radsec from nipype.workflows.dmri.fsl.utils)
-    # pre_fugue = pe.Node(fsl.FUGUE(save_fmap=True),
-    #                     name='ComputeFieldmapFUGUE')
+    # pre_fugue = pe.Node(fsl.FUGUE(save_fmap=True), name='ComputeFieldmapFUGUE')
     # rsec2hz (divide by 2pi)
 
     if phasetype == "phasediff":
@@ -124,7 +120,7 @@ def init_phdiff_wf(omp_nthreads, phasetype="phasediff", name='phdiff_wf'):
     elif phasetype == "phase":
         workflow.__desc__ += """\
 The phase difference used for unwarping was calculated using a method developed
-by Mark Elliot (University of Pennsylvania).
+by Mark Elliott [@pncprocessing].
     """
         # Special case for phase1, phase2 images
         meta = pe.MapNode(ReadSidecarJSON(), name='meta', mem_gb=0.01,

From 2c52d4973c7e31f27f84bfc91795b31f8df3c260 Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Tue, 20 Nov 2018 16:19:33 -0500
Subject: [PATCH 23/34] pylint error

---
 fmriprep/interfaces/fmap.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/fmriprep/interfaces/fmap.py b/fmriprep/interfaces/fmap.py
index ead71a875..f2363212e 100644
--- a/fmriprep/interfaces/fmap.py
+++ b/fmriprep/interfaces/fmap.py
@@ -72,7 +72,7 @@ def _run_interface(self, runtime):
                 mask = sim.binary_erosion(
                     mask, struc,
                     iterations=self.inputs.mask_erode
-                ).astype(np.uint8) # pylint: disable=no-member
+                ).astype(np.uint8)  # pylint: disable=no-member
 
         self._results['out_file'] = fname_presuffix(
             self.inputs.in_file, suffix='_enh', newpath=runtime.cwd)

From 27c51ae5c5f004617ea0b2658425e571208e49db Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Tue, 20 Nov 2018 16:23:38 -0500
Subject: [PATCH 24/34] Fix zenodo

---
 .zenodo.json | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/.zenodo.json b/.zenodo.json
index 5c505d45b..a17fa0a1e 100644
--- a/.zenodo.json
+++ b/.zenodo.json
@@ -88,11 +88,6 @@
       "affiliation": "Department of Psychology, University of Oslo",
       "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",
@@ -103,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",

From f57c0634cacaaf6e81985a2f0bfd377df0d1cc80 Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Tue, 20 Nov 2018 16:52:44 -0500
Subject: [PATCH 25/34] missed a couple formatting changes

---
 fmriprep/interfaces/fmap.py           | 6 ++++--
 fmriprep/workflows/fieldmap/phdiff.py | 3 +--
 2 files changed, 5 insertions(+), 4 deletions(-)

diff --git a/fmriprep/interfaces/fmap.py b/fmriprep/interfaces/fmap.py
index f2363212e..7fc114ab6 100644
--- a/fmriprep/interfaces/fmap.py
+++ b/fmriprep/interfaces/fmap.py
@@ -93,7 +93,8 @@ 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()
 
@@ -200,7 +201,8 @@ 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
 
 
diff --git a/fmriprep/workflows/fieldmap/phdiff.py b/fmriprep/workflows/fieldmap/phdiff.py
index 6b5cc43ca..eb78f4132 100644
--- a/fmriprep/workflows/fieldmap/phdiff.py
+++ b/fmriprep/workflows/fieldmap/phdiff.py
@@ -35,8 +35,7 @@
 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)`
+    magnitude images corresponding to two or more :abbr:`GRE (Gradient Echo sequence)`
     acquisitions. The `original code was taken from nipype
     <https://github.com/nipy/nipype/blob/master/nipype/workflows/dmri/fsl/artifacts.py#L514>`_.
 

From 2f6a294fdc40a07ca7825a47883e4364a23b1db0 Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Sat, 24 Nov 2018 13:25:37 -0500
Subject: [PATCH 26/34] style changes

---
 fmriprep/interfaces/fmap.py | 13 +++++++------
 1 file changed, 7 insertions(+), 6 deletions(-)

diff --git a/fmriprep/interfaces/fmap.py b/fmriprep/interfaces/fmap.py
index 7fc114ab6..64c1a2225 100644
--- a/fmriprep/interfaces/fmap.py
+++ b/fmriprep/interfaces/fmap.py
@@ -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()
 
@@ -201,7 +201,8 @@ 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
 
@@ -496,7 +497,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

From 0f7cdecaa2264721dd2b075715367a5c8c92757d Mon Sep 17 00:00:00 2001
From: Chris Markiewicz <markiewicz@stanford.edu>
Date: Sat, 15 Dec 2018 16:00:31 -0500
Subject: [PATCH 27/34] FIX: Typo in .zenodo.json

---
 .zenodo.json | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.zenodo.json b/.zenodo.json
index 0bd21e030..b2a9d8767 100644
--- a/.zenodo.json
+++ b/.zenodo.json
@@ -107,7 +107,7 @@
       "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",

From f2e861d61782d54bd31331bb701af7e029b6fbbf Mon Sep 17 00:00:00 2001
From: "Christopher J. Markiewicz" <markiewicz@stanford.edu>
Date: Sat, 15 Dec 2018 16:04:41 -0500
Subject: [PATCH 28/34] STY: Revert unneeded style changes

---
 fmriprep/interfaces/fmap.py         | 7 ++++---
 fmriprep/workflows/fieldmap/base.py | 4 +++-
 2 files changed, 7 insertions(+), 4 deletions(-)

diff --git a/fmriprep/interfaces/fmap.py b/fmriprep/interfaces/fmap.py
index 64c1a2225..5a0e4bc05 100644
--- a/fmriprep/interfaces/fmap.py
+++ b/fmriprep/interfaces/fmap.py
@@ -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
 
@@ -260,7 +260,8 @@ 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
 
diff --git a/fmriprep/workflows/fieldmap/base.py b/fmriprep/workflows/fieldmap/base.py
index e6f79d666..93732d99b 100644
--- a/fmriprep/workflows/fieldmap/base.py
+++ b/fmriprep/workflows/fieldmap/base.py
@@ -217,7 +217,9 @@ def init_sdc_wf(fmaps, bold_meta, omp_nthreads=1,
             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,

From d94346f0cbdb5f7d85f0020b018474460cabc665 Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Tue, 18 Dec 2018 12:25:03 -0500
Subject: [PATCH 29/34] Interpret phase image values directly

---
 fmriprep/interfaces/fmap.py           | 10 +++++++---
 fmriprep/workflows/fieldmap/phdiff.py |  4 ++--
 2 files changed, 9 insertions(+), 5 deletions(-)

diff --git a/fmriprep/interfaces/fmap.py b/fmriprep/interfaces/fmap.py
index 5a0e4bc05..29d9d2dd8 100644
--- a/fmriprep/interfaces/fmap.py
+++ b/fmriprep/interfaces/fmap.py
@@ -228,7 +228,8 @@ class Phases2Fieldmap(SimpleInterface):
 
     def _run_interface(self, runtime):
         # Get the echo times
-        fmap_file, merged_metadata = phases2fmap(self.inputs.phase_files, self.inputs.metadatas)
+        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
@@ -549,10 +550,13 @@ def phases2fmap(phase_files, metadatas, newpath=None):
     phase1 = image1.get_fdata()
 
     def rescale_image(img):
+        mask = img > 0
         imax = img.max()
         imin = img.min()
-        scaled = 2 * ((img - imin) / (imax - imin) - 0.5)
-        return np.pi * scaled
+        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)
diff --git a/fmriprep/workflows/fieldmap/phdiff.py b/fmriprep/workflows/fieldmap/phdiff.py
index b25deb027..4b8ae975d 100644
--- a/fmriprep/workflows/fieldmap/phdiff.py
+++ b/fmriprep/workflows/fieldmap/phdiff.py
@@ -117,8 +117,8 @@ def init_phdiff_wf(omp_nthreads, phasetype="phasediff", name='phdiff_wf'):
 
     elif phasetype == "phase":
         workflow.__desc__ += """\
-The phase difference used for unwarping was calculated using a method developed
-by Mark Elliott [@pncprocessing].
+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,

From 6c7ec33664e2c3f8341f4b109c04109f35887196 Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Thu, 20 Dec 2018 12:27:21 -0500
Subject: [PATCH 30/34] Check that encoding is not specified as Bipolar

---
 fmriprep/workflows/fieldmap/base.py | 15 +++++++++++++++
 1 file changed, 15 insertions(+)

diff --git a/fmriprep/workflows/fieldmap/base.py b/fmriprep/workflows/fieldmap/base.py
index 6e4363874..79cd5879b 100644
--- a/fmriprep/workflows/fieldmap/base.py
+++ b/fmriprep/workflows/fieldmap/base.py
@@ -215,6 +215,21 @@ def init_sdc_wf(fmaps, bold_meta, omp_nthreads=1,
             if fmap['type'] == 'phasediff':
                 fmap_estimator_wf.inputs.inputnode.phasediff = fmap['phasediff']
             elif fmap['type'] == '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())

From c27d83931092758d20adafaa185470cb67e90e40 Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Sat, 11 May 2019 15:19:52 -0400
Subject: [PATCH 31/34] change mask name to ds_report_fmap_mask

---
 fmriprep/workflows/fieldmap/phdiff.py | 22 +++++++++++++++-------
 1 file changed, 15 insertions(+), 7 deletions(-)

diff --git a/fmriprep/workflows/fieldmap/phdiff.py b/fmriprep/workflows/fieldmap/phdiff.py
index 4b8ae975d..9212253b0 100644
--- a/fmriprep/workflows/fieldmap/phdiff.py
+++ b/fmriprep/workflows/fieldmap/phdiff.py
@@ -15,7 +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
+8.9.2 in BIDS 1.0.0: two phases and at least one magnitude image
 
 """
 
@@ -31,7 +31,7 @@
 from ...interfaces import Phasediff2Fieldmap, Phases2Fieldmap, DerivativesDataSink
 
 
-def init_phdiff_wf(omp_nthreads, phasetype="phasediff", 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)`
@@ -78,8 +78,9 @@ def init_phdiff_wf(omp_nthreads, phasetype="phasediff", name='phdiff_wf'):
                  name='n4', n_procs=omp_nthreads)
     bet = pe.Node(BETRPT(generate_report=True, frac=0.6, mask=True),
                   name='bet')
-    ds_fmap_mask = pe.Node(DerivativesDataSink(suffix='fmap_mask'), name='ds_report_fmap_mask',
-                           mem_gb=0.01, run_without_submitting=True)
+    ds_report_fmap_mask = pe.Node(DerivativesDataSink(
+        desc='brain', suffix='mask'), name='ds_report_fmap_mask',
+        mem_gb=0.01, run_without_submitting=True)
     # uses mask from bet; outputs a mask
     # dilate = pe.Node(fsl.maths.MathsCommand(
     #     nan2zeros=True, args='-kernel sphere 5 -dilM'), name='MskDilate')
@@ -102,6 +103,9 @@ def init_phdiff_wf(omp_nthreads, phasetype="phasediff", name='phdiff_wf'):
     # 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')
@@ -112,7 +116,7 @@ def init_phdiff_wf(omp_nthreads, phasetype="phasediff", name='phdiff_wf'):
             (meta, compfmap, [('out_dict', 'metadata')]),
             (inputnode, pha2rads, [('phasediff', 'in_file')]),
             (pha2rads, prelude, [('out', 'phase_file')]),
-            (inputnode, ds_fmap_mask, [('phasediff', 'source_file')]),
+            (inputnode, ds_report_fmap_mask, [('phasediff', 'source_file')]),
         ])
 
     elif phasetype == "phase":
@@ -129,7 +133,7 @@ def init_phdiff_wf(omp_nthreads, phasetype="phasediff", name='phdiff_wf'):
             (inputnode, phases2fmap, [('phasediff', 'phase_files')]),
             (phases2fmap, prelude, [('out_file', 'phase_file')]),
             (phases2fmap, compfmap, [('phasediff_metadata', 'metadata')]),
-            (phases2fmap, ds_fmap_mask, [('out_file', 'source_file')])
+            (phases2fmap, ds_report_fmap_mask, [('out_file', 'source_file')])
         ])
 
     workflow.connect([
@@ -139,6 +143,9 @@ def init_phdiff_wf(omp_nthreads, phasetype="phasediff", name='phdiff_wf'):
         (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')]),
@@ -147,7 +154,8 @@ def init_phdiff_wf(omp_nthreads, phasetype="phasediff", name='phdiff_wf'):
         (compfmap, outputnode, [('out_file', 'fmap')]),
         (bet, outputnode, [('mask_file', 'fmap_mask'),
                            ('out_file', 'fmap_ref')]),
-        (bet, ds_fmap_mask, [('out_report', 'in_file')]),
+        (inputnode, ds_report_fmap_mask, [('phasediff', 'source_file')]),
+        (bet, ds_report_fmap_mask, [('out_report', 'in_file')]),
     ])
 
     return workflow

From 4d371a117a07dc36ad7286835914a5e3dd204ddc Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Sat, 11 May 2019 21:34:30 -0400
Subject: [PATCH 32/34] replate type with suffix

---
 fmriprep/interfaces/fmap.py           |  2 +-
 fmriprep/workflows/fieldmap/base.py   | 12 ++++++------
 fmriprep/workflows/fieldmap/phdiff.py |  4 ----
 3 files changed, 7 insertions(+), 11 deletions(-)

diff --git a/fmriprep/interfaces/fmap.py b/fmriprep/interfaces/fmap.py
index 29d9d2dd8..0b128fe27 100644
--- a/fmriprep/interfaces/fmap.py
+++ b/fmriprep/interfaces/fmap.py
@@ -580,7 +580,7 @@ def rescale_image(img):
 
 
 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.
diff --git a/fmriprep/workflows/fieldmap/base.py b/fmriprep/workflows/fieldmap/base.py
index 08d8d0628..b60449cf2 100644
--- a/fmriprep/workflows/fieldmap/base.py
+++ b/fmriprep/workflows/fieldmap/base.py
@@ -200,8 +200,8 @@ def init_sdc_wf(fmaps, bold_meta, omp_nthreads=1,
         ])
 
     # FIELDMAP path
-    if fmap['type'] in ['fieldmap', 'phasediff', 'phase']:
-        outputnode.inputs.method = 'FMB (%s-based)' % fmap['type']
+    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.
         if fmap['suffix'] == 'fieldmap':
@@ -213,14 +213,14 @@ 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'] in ('phasediff', 'phase'):
+        if fmap['suffix'] in ('phasediff', 'phase'):
             from .phdiff import init_phdiff_wf
             fmap_estimator_wf = init_phdiff_wf(omp_nthreads=omp_nthreads,
-                                               phasetype=fmap['type'])
+                                               phasetype=fmap['suffix'])
             # set inputs
-            if fmap['type'] == 'phasediff':
+            if fmap['suffix'] == 'phasediff':
                 fmap_estimator_wf.inputs.inputnode.phasediff = fmap['phasediff']
-            elif fmap['type'] == 'phase':
+            elif fmap['suffix'] == 'phase':
                 # Check that fieldmap is not bipolar
                 fmap_polarity = fmap['metadata'].get('DiffusionScheme', None)
                 if fmap_polarity == 'Bipolar':
diff --git a/fmriprep/workflows/fieldmap/phdiff.py b/fmriprep/workflows/fieldmap/phdiff.py
index 9212253b0..c233f4059 100644
--- a/fmriprep/workflows/fieldmap/phdiff.py
+++ b/fmriprep/workflows/fieldmap/phdiff.py
@@ -143,9 +143,6 @@ def init_phdiff_wf(omp_nthreads, phasetype='phasediff', name='phdiff_wf'):
         (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')]),
@@ -154,7 +151,6 @@ def init_phdiff_wf(omp_nthreads, phasetype='phasediff', name='phdiff_wf'):
         (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')]),
     ])
 

From 9c5b46bf8236d2896bc9419ffb9d92b13164d4fc Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Sat, 11 May 2019 22:09:44 -0400
Subject: [PATCH 33/34] more changes of type to suffix

---
 fmriprep/workflows/bold/base.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/fmriprep/workflows/bold/base.py b/fmriprep/workflows/bold/base.py
index 7d90f9b12..6ce2d9988 100755
--- a/fmriprep/workflows/bold/base.py
+++ b/fmriprep/workflows/bold/base.py
@@ -336,10 +336,10 @@ def init_func_preproc_wf(
         if 'fieldmaps' not in ignore:
             fmaps = layout.get_fieldmap(ref_file, return_list=True)
             for fmap in fmaps:
-                if fmap['type'] == 'phase':
+                if fmap['suffix'] == 'phase':
                     fmap_key = 'phase1'
                 else:
-                    fmap_key = fmap['type']
+                    fmap_key = fmap['suffix']
                 fmap['metadata'] = layout.get_metadata(fmap[fmap_key])
 
         # Run SyN if forced or in the absence of fieldmap correction

From 85fcb4e093ddf26386839d1d4d04c2e9ef8af921 Mon Sep 17 00:00:00 2001
From: Matt Cieslak <mattcieslak@gmail.com>
Date: Mon, 13 May 2019 14:31:41 -0400
Subject: [PATCH 34/34] Support negative numbers in phase images (happens at
 7T)

---
 fmriprep/interfaces/fmap.py | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/fmriprep/interfaces/fmap.py b/fmriprep/interfaces/fmap.py
index 0b128fe27..71c96cc25 100644
--- a/fmriprep/interfaces/fmap.py
+++ b/fmriprep/interfaces/fmap.py
@@ -550,6 +550,13 @@ def phases2fmap(phase_files, metadatas, newpath=None):
     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()