fix: add recenter node, improve documentation

Author: oesteban

sdcflows/interfaces/fmap.py

@@ -187,7 +187,16 @@ class _Phasediff2FieldmapOutputSpec(TraitedSpec):
 
 
 class Phasediff2Fieldmap(SimpleInterface):
-    """Convert a phase difference map into a fieldmap in Hz."""
+    """
+    Convert a phase difference map into a fieldmap in Hz.
+
+    This interface is equivalent to running the following steps:
+      #. Convert from rad to rad/s
+         (``niflow.nipype1.workflows.dmri.fsl.utils.rads2radsec``)
+      #. FUGUE execution: fsl.FUGUE(save_fmap=True)
+      #. Conversion from rad/s to Hz (divide by 2pi, ``rsec2hz``).
+
+    """
 
     input_spec = _Phasediff2FieldmapInputSpec
     output_spec = _Phasediff2FieldmapOutputSpec

sdcflows/workflows/phdiff.py

@@ -29,12 +29,22 @@
 
 
 def init_phdiff_wf(omp_nthreads, name='phdiff_wf'):
-    """
+    r"""
     Distortion correction of EPI sequences using phase-difference maps.
 
     Estimates the fieldmap using a phase-difference image and one or more
     magnitude images corresponding to two or more :abbr:`GRE (Gradient Echo sequence)`
-    acquisitions. The `original code was taken from nipype
+    acquisitions.
+    The most delicate bit of this workflow is the phase-unwrapping process: phase maps
+    are clipped in the range :math:`[0 \dotsb 2 \cdot \pi )`.
+    To find the integer number of offsets that make a region continously smooth with
+    its neighbour, FSL PRELUDE is run [Jenkinson1998]_.
+    FSL PRELUDE takes wrapped maps in the range 0 to 6.28, `as per the user guide
+    <https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FUGUE/Guide#Step_2_-_Getting_.28wrapped.29_phase_in_radians>`__.
+    For the phase-difference maps, recentering back to :math:`[-\pi \dotsb \pi )` is necessary.
+    After some massaging and the application of the effective echo spacing parameter,
+    the phase-difference maps can be converted into a *B0 field map* in Hz units.
+    The `original code was taken from nipype
     <https://github.com/nipy/nipype/blob/0.12.1/nipype/workflows/dmri/fsl/artifacts.py#L514>`_.
 
     Workflow Graph
@@ -98,11 +108,15 @@ def init_phdiff_wf(omp_nthreads, name='phdiff_wf'):
     #     nan2zeros=True, args='-kernel sphere 5 -dilM'), name='MskDilate')
 
     # phase diff -> radians
-    phmap2rads = pe.Node(PhaseMap2rads(), name='phmap2rads')
+    phmap2rads = pe.Node(PhaseMap2rads(), name='phmap2rads',
+                         run_without_submitting=True)
 
     # FSL PRELUDE will perform phase-unwrapping
     prelude = pe.Node(fsl.PRELUDE(), name='prelude')
 
+    recenter = pe.Node(niu.Function(function=_recenter),
+                       name='recenter', run_without_submitting=True)
+
     denoise = pe.Node(fsl.SpatialFilter(operation='median', kernel_shape='sphere',
                                         kernel_size=3), name='denoise')
 
@@ -112,11 +126,6 @@ def init_phdiff_wf(omp_nthreads, name='phdiff_wf'):
 
     compfmap = pe.Node(Phasediff2Fieldmap(), name='compfmap')
 
-    # The phdiff2fmap interface is equivalent to:
-    # rad2rsec (using rads2radsec from niflow.nipype1.workflows.dmri.fsl.utils)
-    # pre_fugue = pe.Node(fsl.FUGUE(save_fmap=True), name='ComputeFieldmapFUGUE')
-    # rsec2hz (divide by 2pi)
-
     workflow.connect([
         (inputnode, compfmap, [('metadata', 'metadata')]),
         (inputnode, magmrg, [('magnitude', 'in_files')]),
@@ -126,7 +135,8 @@ def init_phdiff_wf(omp_nthreads, name='phdiff_wf'):
         (bet, prelude, [('mask_file', 'mask_file')]),
         (inputnode, phmap2rads, [('phasediff', 'in_file')]),
         (phmap2rads, prelude, [('out_file', 'phase_file')]),
-        (prelude, denoise, [('unwrapped_phase_file', 'in_file')]),
+        (prelude, recenter, [('unwrapped_phase_file', 'in_file')]),
+        (recenter, denoise, [('out', 'in_file')]),
         (denoise, demean, [('out_file', 'in_file')]),
         (demean, cleanup_wf, [('out', 'inputnode.in_file')]),
         (bet, cleanup_wf, [('mask_file', 'inputnode.in_mask')]),
@@ -137,3 +147,21 @@ def init_phdiff_wf(omp_nthreads, name='phdiff_wf'):
     ])
 
     return workflow
+
+
+def _recenter(in_file, offset=None):
+    """Recenter the phase-map distribution to the -pi..pi range."""
+    from os.path import basename
+    import numpy as np
+    import nibabel as nb
+    from nipype.utils.filemanip import fname_presuffix
+
+    if offset is None:
+        offset = np.pi
+
+    nii = nb.load(in_file)
+    data = nii.get_fdata(dtype='float32') - offset
+
+    out_file = fname_presuffix(basename(in_file), suffix='_recentered')
+    nb.Nifti1Image(data, nii.affine, nii.header).to_filename(out_file)
+    return out_file