ENH: Stop using siemens2rads from old nipype workflows

A new interface is proposed, to be used in phase-difference based
workflows.

@Aksoo identified a potential issue with siemens2rads
(#30 (comment)),
although the re-centering did not seem to make a big difference.

This implementation of the rescaling uses percentiles to robustly
calculate the range of the input image (as recommended in the original
paper about FSL PRELUDE), and makes sure the output data type is
adequate (float32).

Author: oesteban

sdcflows/interfaces/fmap.py

@@ -200,6 +200,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 +552,22 @@ 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."""
+    im = nb.load(in_file)
+    data = im.get_fdata().astype('float32')
+    hdr = im.header.copy()
+
+    data -= np.percentile(data, 2)
+    data[data < 0] = 0
+    data = 2.0 * np.pi * data / np.percentile(data, 98)
+    data[data > 2.0 * np.pi] = 2.0 * np.pi
+    hdr.set_data_dtype(np.float32)
+    hdr.set_xyzt_units('mm')
+    hdr['datatype'] = 16
+
+    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,12 +20,12 @@
 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'):
@@ -98,7 +98,7 @@ 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')
 
     # FSL PRELUDE will perform phase-unwrapping
     prelude = pe.Node(fsl.PRELUDE(), name='prelude')
@@ -124,8 +124,8 @@ def init_phdiff_wf(omp_nthreads, 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')]),
+        (inputnode, phmap2rads, [('phasediff', 'in_file')]),
+        (phmap2rads, prelude, [('out_file', 'phase_file')]),
         (prelude, denoise, [('unwrapped_phase_file', 'in_file')]),
         (denoise, demean, [('out_file', 'in_file')]),
         (demean, cleanup_wf, [('out', 'inputnode.in_file')]),