ENH: Stop using siemens2rads from old nipype workflows

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
@@ -200,6 +209,27 @@ def _run_interface(self, runtime):
         return runtime
 
 
+class _PhaseMap2radsInputSpec(BaseInterfaceInputSpec):
+    in_file = File(exists=True, mandatory=True, desc='input (wrapped) phase map')
+
+
+class _PhaseMap2radsOutputSpec(TraitedSpec):
+    out_file = File(desc='the phase map in the range 0 - 6.28')
+
+
+class PhaseMap2rads(SimpleInterface):
+    """Convert a phase map in a.u. to radians."""
+
+    input_spec = _PhaseMap2radsInputSpec
+    output_spec = _PhaseMap2radsOutputSpec
+
+    def _run_interface(self, runtime):
+        self._results['out_file'] = au2rads(
+            self.inputs.in_file,
+            newpath=runtime.cwd)
+        return runtime
+
+
 def _despike2d(data, thres, neigh=None):
     """Despike axial slices, as done in FSL's ``epiunwarp``."""
     if neigh is None:
@@ -531,3 +561,32 @@ 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 au2rads(in_file, newpath=None):
+    """Convert the input phase difference map in arbitrary units (a.u.) to rads."""
+    from scipy.stats import mode
+    im = nb.load(in_file)
+    data = im.get_fdata(dtype='float32')
+    hdr = im.header.copy()
+
+    # First center data around 0.0.
+    data -= mode(data, axis=None)[0][0]
+
+    # Scale lower tail
+    data[data < 0] = - np.pi * data[data < 0] / data[data < 0].min()
+
+    # Scale upper tail
+    data[data > 0] = np.pi * data[data > 0] / data[data > 0].max()
+
+    # Offset to 0 - 2pi
+    data += np.pi
+
+    # Clip
+    data = np.clip(data, 0.0, 2 * np.pi)
+
+    hdr.set_data_dtype(np.float32)
+    hdr.set_xyzt_units('mm')
+    out_file = fname_presuffix(in_file, suffix='_rads', newpath=newpath)
+    nb.Nifti1Image(data, im.affine, hdr).to_filename(out_file)
+    return out_file

sdcflows/workflows/phdiff.py

@@ -20,21 +20,31 @@
 from nipype.interfaces import ants, fsl, utility as niu
 from nipype.pipeline import engine as pe
 from niflow.nipype1.workflows.dmri.fsl.utils import (
-    siemens2rads, demean_image, cleanup_edge_pipeline)
+    demean_image, cleanup_edge_pipeline)
 from niworkflows.engine.workflows import LiterateWorkflow as Workflow
 from niworkflows.interfaces.images import IntraModalMerge
 from niworkflows.interfaces.masks import BETRPT
 
-from ..interfaces.fmap import Phasediff2Fieldmap
+from ..interfaces.fmap import Phasediff2Fieldmap, PhaseMap2rads
 
 
 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 [Jenkinson2003]_.
+    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
@@ -68,6 +78,11 @@ def init_phdiff_wf(omp_nthreads, name='phdiff_wf'):
     fmap : pathlike
         The estimated fieldmap in Hz
 
+    References
+    ----------
+    .. [Jenkinson2003] Jenkinson, M. (2003) Fast, automated, N-dimensional phase-unwrapping
+        algorithm. MRM 49(1):193-197. doi:`10.1002/mrm.10354 <10.1002/mrm.10354>`__.
+
     """
     workflow = Workflow(name=name)
     workflow.__desc__ = """\
@@ -98,11 +113,15 @@ def init_phdiff_wf(omp_nthreads, name='phdiff_wf'):
     #     nan2zeros=True, args='-kernel sphere 5 -dilM'), name='MskDilate')
 
     # phase diff -> radians
-    pha2rads = pe.Node(niu.Function(function=siemens2rads), name='pha2rads')
+    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,21 +131,17 @@ 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')]),
         (magmrg, n4, [('out_avg', 'input_image')]),
         (n4, prelude, [('output_image', 'magnitude_file')]),
         (n4, bet, [('output_image', 'in_file')]),
         (bet, prelude, [('mask_file', 'mask_file')]),
-        (inputnode, pha2rads, [('phasediff', 'in_file')]),
-        (pha2rads, prelude, [('out', 'phase_file')]),
-        (prelude, denoise, [('unwrapped_phase_file', 'in_file')]),
+        (inputnode, phmap2rads, [('phasediff', 'in_file')]),
+        (phmap2rads, prelude, [('out_file', 'phase_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 +152,22 @@ def init_phdiff_wf(omp_nthreads, name='phdiff_wf'):
     ])
 
     return workflow
+
+
+def _recenter(in_file):
+    """Recenter the phase-map distribution to the -pi..pi range."""
+    from os import getcwd
+    import numpy as np
+    import nibabel as nb
+    from nipype.utils.filemanip import fname_presuffix
+
+    nii = nb.load(in_file)
+    data = nii.get_fdata(dtype='float32')
+    msk = data != 0
+    msk[data == 0] = False
+    data[msk] -= np.median(data[msk])
+
+    out_file = fname_presuffix(in_file, suffix='_recentered',
+                               newpath=getcwd())
+    nb.Nifti1Image(data, nii.affine, nii.header).to_filename(out_file)
+    return out_file