fix: implement a bilateral scaling around the mode, better recentering

Author: oesteban

sdcflows/interfaces/fmap.py

@@ -565,17 +565,28 @@ def _delta_te(in_values, te1=None, te2=None):
 
 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()
 
-    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
+    # First center data around 0.0.
+    data -= mode(data, axis=None)[0][0]
+
+    # Scale lower tail
+    data[data < 0] = - np.pi * data[data < 0] / np.percentile(data[data < 0], 2)
+
+    # Scale upper tail
+    data[data > 0] = np.pi * data[data > 0] / np.percentile(data[data > 0], 98)
+
+    # 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

@@ -141,7 +141,6 @@ def init_phdiff_wf(omp_nthreads, name='phdiff_wf'):
         (inputnode, phmap2rads, [('phasediff', 'in_file')]),
         (phmap2rads, prelude, [('out_file', 'phase_file')]),
         (prelude, recenter, [('unwrapped_phase_file', 'in_file')]),
-        (bet, recenter, [('mask_file', 'in_mask')]),
         (recenter, denoise, [('out', 'in_file')]),
         (denoise, demean, [('out_file', 'in_file')]),
         (demean, cleanup_wf, [('out', 'inputnode.in_file')]),
@@ -155,24 +154,18 @@ def init_phdiff_wf(omp_nthreads, name='phdiff_wf'):
     return workflow
 
 
-def _recenter(in_file, in_mask=None, offset=None):
+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
 
-    if offset is None:
-        offset = np.pi
-
     nii = nb.load(in_file)
     data = nii.get_fdata(dtype='float32')
-
-    msk = data > 1.e-6
-    if in_mask is not None:
-        msk = nb.load(in_mask).get_fdata(dtype='float32') > 1.e-4
-
-    data[msk] -= offset
+    msk = data != 0
+    msk[data == 0] = False
+    data[msk] -= np.median(data[msk])
 
     out_file = fname_presuffix(in_file, suffix='_recentered',
                                newpath=getcwd())