From 04637a5a91974372c478eff5cd8644859f65f2a2 Mon Sep 17 00:00:00 2001 From: oesteban Date: Mon, 18 Nov 2019 10:00:12 -0800 Subject: [PATCH 1/8] 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 (https://github.com/poldracklab/sdcflows/pull/30#discussion_r346640162), 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). --- sdcflows/interfaces/fmap.py | 40 ++++++++++++++++++++++++++++++++++++ sdcflows/workflows/phdiff.py | 10 ++++----- 2 files changed, 45 insertions(+), 5 deletions(-) diff --git a/sdcflows/interfaces/fmap.py b/sdcflows/interfaces/fmap.py index 1d9ddff521..1b4b0035ee 100644 --- a/sdcflows/interfaces/fmap.py +++ b/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 diff --git a/sdcflows/workflows/phdiff.py b/sdcflows/workflows/phdiff.py index 676f0718c4..765eebc8ce 100644 --- a/sdcflows/workflows/phdiff.py +++ b/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')]), From e21ed8bba28aa2ade3f6d9ffbf62c742b9687846 Mon Sep 17 00:00:00 2001 From: Oscar Esteban Date: Mon, 18 Nov 2019 13:33:42 -0500 Subject: [PATCH 2/8] Update sdcflows/interfaces/fmap.py [skip ci] Co-Authored-By: Chris Markiewicz --- sdcflows/interfaces/fmap.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sdcflows/interfaces/fmap.py b/sdcflows/interfaces/fmap.py index 1b4b0035ee..37617379af 100644 --- a/sdcflows/interfaces/fmap.py +++ b/sdcflows/interfaces/fmap.py @@ -557,7 +557,7 @@ 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.""" im = nb.load(in_file) - data = im.get_fdata().astype('float32') + data = im.get_fdata(dtype='float32') hdr = im.header.copy() data -= np.percentile(data, 2) From aea890bdb7787160929395c1464a5cf2b3a5a048 Mon Sep 17 00:00:00 2001 From: oesteban Date: Mon, 18 Nov 2019 10:38:59 -0800 Subject: [PATCH 3/8] enh: remove superfluous data type setting --- sdcflows/interfaces/fmap.py | 1 - 1 file changed, 1 deletion(-) diff --git a/sdcflows/interfaces/fmap.py b/sdcflows/interfaces/fmap.py index 37617379af..2a93918961 100644 --- a/sdcflows/interfaces/fmap.py +++ b/sdcflows/interfaces/fmap.py @@ -566,7 +566,6 @@ def au2rads(in_file, newpath=None): 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) From 616c16ffab36d6a6d99be423c205914946ad2db2 Mon Sep 17 00:00:00 2001 From: oesteban Date: Mon, 18 Nov 2019 11:50:48 -0800 Subject: [PATCH 4/8] fix: add recenter node, improve documentation --- sdcflows/interfaces/fmap.py | 11 ++++++++- sdcflows/workflows/phdiff.py | 46 +++++++++++++++++++++++++++++------- 2 files changed, 47 insertions(+), 10 deletions(-) diff --git a/sdcflows/interfaces/fmap.py b/sdcflows/interfaces/fmap.py index 2a93918961..ae43efdc7e 100644 --- a/sdcflows/interfaces/fmap.py +++ b/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 diff --git a/sdcflows/workflows/phdiff.py b/sdcflows/workflows/phdiff.py index 765eebc8ce..ddb7a5cf7e 100644 --- a/sdcflows/workflows/phdiff.py +++ b/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 + `__. + 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 `_. 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 From 97ed3db204d51756a266833c7e656f5c7d80b9ad Mon Sep 17 00:00:00 2001 From: oesteban Date: Mon, 18 Nov 2019 12:20:44 -0800 Subject: [PATCH 5/8] docs: fix missing reference --- sdcflows/workflows/phdiff.py | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/sdcflows/workflows/phdiff.py b/sdcflows/workflows/phdiff.py index ddb7a5cf7e..bd062bb803 100644 --- a/sdcflows/workflows/phdiff.py +++ b/sdcflows/workflows/phdiff.py @@ -38,7 +38,7 @@ def init_phdiff_wf(omp_nthreads, name='phdiff_wf'): 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]_. + its neighbour, FSL PRELUDE is run [Jenkinson2003]_. FSL PRELUDE takes wrapped maps in the range 0 to 6.28, `as per the user guide `__. For the phase-difference maps, recentering back to :math:`[-\pi \dotsb \pi )` is necessary. @@ -78,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__ = """\ From 588d3fdbb880e9fd7a5b185aed503fa2c86769e5 Mon Sep 17 00:00:00 2001 From: oesteban Date: Mon, 18 Nov 2019 12:42:05 -0800 Subject: [PATCH 6/8] fix: apply recenter only within the mask FSL PRELUDE will return some values unwrapped below zero, making it hard to differenciate these pixels from the mask. This commit addresses that particular issue. --- sdcflows/workflows/phdiff.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/sdcflows/workflows/phdiff.py b/sdcflows/workflows/phdiff.py index bd062bb803..56a8d58cb8 100644 --- a/sdcflows/workflows/phdiff.py +++ b/sdcflows/workflows/phdiff.py @@ -141,6 +141,7 @@ 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')]), @@ -154,9 +155,9 @@ def init_phdiff_wf(omp_nthreads, name='phdiff_wf'): return workflow -def _recenter(in_file, offset=None): +def _recenter(in_file, in_mask=None, offset=None): """Recenter the phase-map distribution to the -pi..pi range.""" - from os.path import basename + from os import getcwd import numpy as np import nibabel as nb from nipype.utils.filemanip import fname_presuffix @@ -165,8 +166,15 @@ def _recenter(in_file, offset=None): offset = np.pi nii = nb.load(in_file) - data = nii.get_fdata(dtype='float32') - offset + data = nii.get_fdata(dtype='float32') - out_file = fname_presuffix(basename(in_file), suffix='_recentered') + 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 + + out_file = fname_presuffix(in_file, suffix='_recentered', + newpath=getcwd()) nb.Nifti1Image(data, nii.affine, nii.header).to_filename(out_file) return out_file From 7ec8e91da02ed7b44d194199d771bd557a31eeea Mon Sep 17 00:00:00 2001 From: oesteban Date: Mon, 18 Nov 2019 16:31:09 -0800 Subject: [PATCH 7/8] fix: implement a bilateral scaling around the mode, better recentering --- sdcflows/interfaces/fmap.py | 21 ++++++++++++++++----- sdcflows/workflows/phdiff.py | 15 ++++----------- 2 files changed, 20 insertions(+), 16 deletions(-) diff --git a/sdcflows/interfaces/fmap.py b/sdcflows/interfaces/fmap.py index ae43efdc7e..4ae5f817c9 100644 --- a/sdcflows/interfaces/fmap.py +++ b/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 diff --git a/sdcflows/workflows/phdiff.py b/sdcflows/workflows/phdiff.py index 56a8d58cb8..5b5ad4f011 100644 --- a/sdcflows/workflows/phdiff.py +++ b/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()) From 82a108a92b203e739951946e547c399decdee199 Mon Sep 17 00:00:00 2001 From: oesteban Date: Mon, 18 Nov 2019 22:12:58 -0800 Subject: [PATCH 8/8] fix: address https://github.com/poldracklab/sdcflows/pull/50#issuecomment-555318605 --- sdcflows/interfaces/fmap.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sdcflows/interfaces/fmap.py b/sdcflows/interfaces/fmap.py index 4ae5f817c9..56ff9f8a9e 100644 --- a/sdcflows/interfaces/fmap.py +++ b/sdcflows/interfaces/fmap.py @@ -574,10 +574,10 @@ def au2rads(in_file, newpath=None): data -= mode(data, axis=None)[0][0] # Scale lower tail - data[data < 0] = - np.pi * data[data < 0] / np.percentile(data[data < 0], 2) + data[data < 0] = - np.pi * data[data < 0] / data[data < 0].min() # Scale upper tail - data[data > 0] = np.pi * data[data > 0] / np.percentile(data[data > 0], 98) + data[data > 0] = np.pi * data[data > 0] / data[data > 0].max() # Offset to 0 - 2pi data += np.pi