diff --git a/docs/api.rst b/docs/api.rst index dbe8a0e76a..89fc44eefe 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -8,5 +8,6 @@ Information on specific functions, classes, and methods. api/sdcflows.fieldmaps api/sdcflows.interfaces + api/sdcflows.utils api/sdcflows.viz api/sdcflows.workflows \ No newline at end of file diff --git a/sdcflows/interfaces/fmap.py b/sdcflows/interfaces/fmap.py index 9f2d1cb7a7..f9e804c89b 100644 --- a/sdcflows/interfaces/fmap.py +++ b/sdcflows/interfaces/fmap.py @@ -13,214 +13,88 @@ """ -import numpy as np -import nibabel as nb from nipype import logging -from nipype.utils.filemanip import fname_presuffix from nipype.interfaces.base import ( - BaseInterfaceInputSpec, TraitedSpec, File, isdefined, traits, - SimpleInterface) + BaseInterfaceInputSpec, + TraitedSpec, + File, + traits, + SimpleInterface, +) -LOGGER = logging.getLogger('nipype.interface') +LOGGER = logging.getLogger("nipype.interface") -class _SubtractPhasesInputSpec(BaseInterfaceInputSpec): - in_phases = traits.List(File(exists=True), min=1, max=2, - desc='input phase maps') - in_meta = traits.List(traits.Dict(), min=1, max=2, - desc='metadata corresponding to the inputs') - - -class _SubtractPhasesOutputSpec(TraitedSpec): - phase_diff = File(exists=True, desc='phase difference map') - metadata = traits.Dict(desc='output metadata') - - -class SubtractPhases(SimpleInterface): - """Calculate a phase difference map.""" - - input_spec = _SubtractPhasesInputSpec - output_spec = _SubtractPhasesOutputSpec - - def _run_interface(self, runtime): - if len(self.inputs.in_phases) != len(self.inputs.in_meta): - raise ValueError( - 'Length of input phase-difference maps and metadata files ' - 'should match.') - - if len(self.inputs.in_phases) == 1: - self._results['phase_diff'] = self.inputs.in_phases[0] - self._results['metadata'] = self.inputs.in_meta[0] - return runtime - - self._results['phase_diff'], self._results['metadata'] = \ - _subtract_phases(self.inputs.in_phases, - self.inputs.in_meta, - newpath=runtime.cwd) - - return runtime - - -class _FieldEnhanceInputSpec(BaseInterfaceInputSpec): - in_file = File(exists=True, mandatory=True, desc='input fieldmap') - in_mask = File(exists=True, desc='brain mask') - in_magnitude = File(exists=True, desc='input magnitude') - unwrap = traits.Bool(False, usedefault=True, desc='run phase unwrap') - despike = traits.Bool(True, usedefault=True, desc='run despike filter') - bspline_smooth = traits.Bool(True, usedefault=True, desc='run 3D bspline smoother') - mask_erode = traits.Int(1, usedefault=True, desc='mask erosion iterations') - despike_threshold = traits.Float(0.2, usedefault=True, desc='mask erosion iterations') - num_threads = traits.Int(1, usedefault=True, nohash=True, desc='number of jobs') +class _PhaseMap2radsInputSpec(BaseInterfaceInputSpec): + in_file = File(exists=True, mandatory=True, desc="input (wrapped) phase map") -class _FieldEnhanceOutputSpec(TraitedSpec): - out_file = File(desc='the output fieldmap') - out_unwrapped = File(desc='unwrapped fieldmap') +class _PhaseMap2radsOutputSpec(TraitedSpec): + out_file = File(desc="the phase map in the range 0 - 6.28") -class FieldEnhance(SimpleInterface): - """Massage the input fieldmap (masking, despiking, etc.).""" +class PhaseMap2rads(SimpleInterface): + """Convert a phase map given in a.u. (e.g., 0-4096) to radians.""" - input_spec = _FieldEnhanceInputSpec - output_spec = _FieldEnhanceOutputSpec + input_spec = _PhaseMap2radsInputSpec + output_spec = _PhaseMap2radsOutputSpec def _run_interface(self, runtime): - from scipy import ndimage as sim - - fmap_nii = nb.load(self.inputs.in_file) - data = np.squeeze(fmap_nii.get_fdata(dtype='float32')) - - # Despike / denoise (no-mask) - if self.inputs.despike: - data = _despike2d(data, self.inputs.despike_threshold) - - mask = None - if isdefined(self.inputs.in_mask): - masknii = nb.load(self.inputs.in_mask) - mask = np.asanyarray(masknii.dataobj).astype('uint8') - - # Dilate mask - if self.inputs.mask_erode > 0: - struc = sim.iterate_structure(sim.generate_binary_structure(3, 2), 1) - mask = sim.binary_erosion( - mask, struc, - iterations=self.inputs.mask_erode - ).astype(np.uint8) # pylint: disable=no-member - - self._results['out_file'] = fname_presuffix( - self.inputs.in_file, suffix='_enh', newpath=runtime.cwd) - datanii = nb.Nifti1Image(data, fmap_nii.affine, fmap_nii.header) - - if self.inputs.unwrap: - data = _unwrap(data, self.inputs.in_magnitude, mask) - self._results['out_unwrapped'] = fname_presuffix( - self.inputs.in_file, suffix='_unwrap', newpath=runtime.cwd) - nb.Nifti1Image(data, fmap_nii.affine, fmap_nii.header).to_filename( - self._results['out_unwrapped']) - - if not self.inputs.bspline_smooth: - datanii.to_filename(self._results['out_file']) - return runtime - else: - from ..utils import bspline as fbsp - from statsmodels.robust.scale import mad - - # Fit BSplines (coarse) - bspobj = fbsp.BSplineFieldmap(datanii, weights=mask, - njobs=self.inputs.num_threads) - bspobj.fit() - smoothed1 = bspobj.get_smoothed() - - # Manipulate the difference map - diffmap = data - smoothed1.get_fdata(dtype='float32') - sderror = mad(diffmap[mask > 0]) - LOGGER.info('SD of error after B-Spline fitting is %f', sderror) - errormask = np.zeros_like(diffmap) - errormask[np.abs(diffmap) > (10 * sderror)] = 1 - errormask *= mask - - nslices = 0 - try: - errorslice = np.squeeze(np.argwhere(errormask.sum(0).sum(0) > 0)) - nslices = errorslice[-1] - errorslice[0] - except IndexError: # mask is empty, do not refine - pass - - if nslices > 1: - diffmapmsk = mask[..., errorslice[0]:errorslice[-1]] - diffmapnii = nb.Nifti1Image( - diffmap[..., errorslice[0]:errorslice[-1]] * diffmapmsk, - datanii.affine, datanii.header) - - bspobj2 = fbsp.BSplineFieldmap(diffmapnii, knots_zooms=[24., 24., 4.], - njobs=self.inputs.num_threads) - bspobj2.fit() - smoothed2 = bspobj2.get_smoothed().get_fdata(dtype='float32') - - final = smoothed1.get_fdata(dtype='float32').copy() - final[..., errorslice[0]:errorslice[-1]] += smoothed2 - else: - final = smoothed1.get_fdata(dtype='float32') - - nb.Nifti1Image(final, datanii.affine, datanii.header).to_filename( - self._results['out_file']) + from ..utils.phasemanip import au2rads + self._results["out_file"] = au2rads(self.inputs.in_file, newpath=runtime.cwd) return runtime -class _FieldToRadSInputSpec(BaseInterfaceInputSpec): - in_file = File(exists=True, mandatory=True, desc='input fieldmap') - fmap_range = traits.Float(desc='range of input field map') +class _SubtractPhasesInputSpec(BaseInterfaceInputSpec): + in_phases = traits.List(File(exists=True), min=1, max=2, desc="input phase maps") + in_meta = traits.List( + traits.Dict(), min=1, max=2, desc="metadata corresponding to the inputs" + ) -class _FieldToRadSOutputSpec(TraitedSpec): - out_file = File(desc='the output fieldmap') - fmap_range = traits.Float(desc='range of input field map') +class _SubtractPhasesOutputSpec(TraitedSpec): + phase_diff = File(exists=True, desc="phase difference map") + metadata = traits.Dict(desc="output metadata") -class FieldToRadS(SimpleInterface): - """Convert from arbitrary units to rad/s.""" +class SubtractPhases(SimpleInterface): + """Calculate a phase difference map.""" - input_spec = _FieldToRadSInputSpec - output_spec = _FieldToRadSOutputSpec + input_spec = _SubtractPhasesInputSpec + output_spec = _SubtractPhasesOutputSpec def _run_interface(self, runtime): - fmap_range = None - if isdefined(self.inputs.fmap_range): - fmap_range = self.inputs.fmap_range - self._results['out_file'], self._results['fmap_range'] = _torads( - self.inputs.in_file, fmap_range, newpath=runtime.cwd) - return runtime - - -class _FieldToHzInputSpec(BaseInterfaceInputSpec): - in_file = File(exists=True, mandatory=True, desc='input fieldmap') - range_hz = traits.Float(mandatory=True, desc='range of input field map') - - -class _FieldToHzOutputSpec(TraitedSpec): - out_file = File(desc='the output fieldmap') + if len(self.inputs.in_phases) != len(self.inputs.in_meta): + raise ValueError( + "Length of input phase-difference maps and metadata files " + "should match." + ) + if len(self.inputs.in_phases) == 1: + self._results["phase_diff"] = self.inputs.in_phases[0] + self._results["metadata"] = self.inputs.in_meta[0] + return runtime -class FieldToHz(SimpleInterface): - """Convert from arbitrary units to Hz.""" + from ..utils.phasemanip import subtract_phases as _subtract_phases - input_spec = _FieldToHzInputSpec - output_spec = _FieldToHzOutputSpec + # Discard in_meta traits with copy(), so that pop() works. + self._results["phase_diff"], self._results["metadata"] = _subtract_phases( + self.inputs.in_phases, + (self.inputs.in_meta[0].copy(), self.inputs.in_meta[1].copy()), + newpath=runtime.cwd, + ) - def _run_interface(self, runtime): - self._results['out_file'] = _tohz( - self.inputs.in_file, self.inputs.range_hz, newpath=runtime.cwd) return runtime class _Phasediff2FieldmapInputSpec(BaseInterfaceInputSpec): - in_file = File(exists=True, mandatory=True, desc='input fieldmap') - metadata = traits.Dict(mandatory=True, desc='BIDS metadata dictionary') + in_file = File(exists=True, mandatory=True, desc="input fieldmap") + metadata = traits.Dict(mandatory=True, desc="BIDS metadata dictionary") class _Phasediff2FieldmapOutputSpec(TraitedSpec): - out_file = File(desc='the output fieldmap') + out_file = File(desc="the output fieldmap") class Phasediff2Fieldmap(SimpleInterface): @@ -239,242 +113,26 @@ class Phasediff2Fieldmap(SimpleInterface): output_spec = _Phasediff2FieldmapOutputSpec def _run_interface(self, runtime): - self._results['out_file'] = phdiff2fmap( - self.inputs.in_file, - _delta_te(self.inputs.metadata), - newpath=runtime.cwd) - 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 - - -class _FUGUEvsm2ANTSwarpInputSpec(BaseInterfaceInputSpec): - in_file = File(exists=True, mandatory=True, - desc='input displacements field map') - pe_dir = traits.Enum('i', 'i-', 'j', 'j-', 'k', 'k-', - desc='phase-encoding axis') - - -class _FUGUEvsm2ANTSwarpOutputSpec(TraitedSpec): - out_file = File(desc='the output warp field') - fieldmap = File(desc='field map in mm') - - -class FUGUEvsm2ANTSwarp(SimpleInterface): - """Convert a voxel-shift-map to ants warp.""" - - _dtype = ' 1e-6 and - (abs(thisval - patch_med) / patch_range) > thres): - data[i, j, k] = patch_med - return data - - -def _unwrap(fmap_data, mag_file, mask=None): - from math import pi - from nipype.interfaces.fsl import PRELUDE - magnii = nb.load(mag_file) - - if mask is None: - mask = np.ones_like(fmap_data, dtype=np.uint8) - - fmapmax = max(abs(fmap_data[mask > 0].min()), fmap_data[mask > 0].max()) - fmap_data *= pi / fmapmax - - nb.Nifti1Image(fmap_data, magnii.affine).to_filename('fmap_rad.nii.gz') - nb.Nifti1Image(mask, magnii.affine).to_filename('fmap_mask.nii.gz') - nb.Nifti1Image(magnii.get_fdata(dtype='float32'), - magnii.affine).to_filename('fmap_mag.nii.gz') - - # Run prelude - res = PRELUDE(phase_file='fmap_rad.nii.gz', - magnitude_file='fmap_mag.nii.gz', - mask_file='fmap_mask.nii.gz').run() - - unwrapped = nb.load( - res.outputs.unwrapped_phase_file).get_fdata(dtype='float32') * (fmapmax / pi) - return unwrapped - - -def _torads(in_file, fmap_range=None, newpath=None): - """ - Convert a field map to rad/s units. - - If fmap_range is None, the range of the fieldmap - will be automatically calculated. - - Use fmap_range=0.5 to convert from Hz to rad/s - - """ - from math import pi - import nibabel as nb - from nipype.utils.filemanip import fname_presuffix - - out_file = fname_presuffix(in_file, suffix='_rad', newpath=newpath) - fmapnii = nb.load(in_file) - fmapdata = fmapnii.get_fdata(dtype='float32') - - if fmap_range is None: - fmap_range = max(abs(fmapdata.min()), fmapdata.max()) - fmapdata = fmapdata * (pi / fmap_range) - out_img = nb.Nifti1Image(fmapdata, fmapnii.affine, fmapnii.header) - out_img.set_data_dtype('float32') - out_img.to_filename(out_file) - return out_file, fmap_range - - -def _tohz(in_file, range_hz, newpath=None): - """Convert a field map to Hz units.""" - from math import pi - import nibabel as nb - from nipype.utils.filemanip import fname_presuffix - - out_file = fname_presuffix(in_file, suffix='_hz', newpath=newpath) - fmapnii = nb.load(in_file) - fmapdata = fmapnii.get_fdata(dtype='float32') - fmapdata = fmapdata * (range_hz / pi) - out_img = nb.Nifti1Image(fmapdata, fmapnii.affine, fmapnii.header) - out_img.set_data_dtype('float32') - out_img.to_filename(out_file) - return out_file - - -def phdiff2fmap(in_file, delta_te, newpath=None): - r""" - Convert the input phase-difference map into a fieldmap in Hz. - - Uses eq. (1) of [Hutton2002]_: - - .. math:: - - \Delta B_0 (\text{T}^{-1}) = \frac{\Delta \Theta}{2\pi\gamma \Delta\text{TE}} - - - In this case, we do not take into account the gyromagnetic ratio of the - proton (:math:`\gamma`), since it will be applied inside TOPUP: - - .. math:: - - \Delta B_0 (\text{Hz}) = \frac{\Delta \Theta}{2\pi \Delta\text{TE}} - - References - ---------- - .. [Hutton2002] Hutton et al., Image Distortion Correction in fMRI: A Quantitative - Evaluation, NeuroImage 16(1):217-240, 2002. doi:`10.1006/nimg.2001.1054 - `_. - - - """ - import math - import numpy as np - import nibabel as nb - from nipype.utils.filemanip import fname_presuffix - # GYROMAG_RATIO_H_PROTON_MHZ = 42.576 - - out_file = fname_presuffix(in_file, suffix='_fmap', newpath=newpath) - image = nb.load(in_file) - data = (image.get_fdata(dtype='float32') / (2. * math.pi * delta_te)) - nii = nb.Nifti1Image(data, image.affine, image.header) - nii.set_data_dtype(np.float32) - nii.to_filename(out_file) - return out_file - - def _delta_te(in_values, te1=None, te2=None): r"""Read :math:`\Delta_\text{TE}` from BIDS metadata dict.""" if isinstance(in_values, float): te2 = in_values - te1 = 0. + te1 = 0.0 if isinstance(in_values, dict): - te1 = in_values.get('EchoTime1') - te2 = in_values.get('EchoTime2') + te1 = in_values.get("EchoTime1") + te2 = in_values.get("EchoTime2") if not all((te1, te2)): - te2 = in_values.get('EchoTimeDifference') + te2 = in_values.get("EchoTimeDifference") te1 = 0 if isinstance(in_values, list): @@ -486,69 +144,17 @@ def _delta_te(in_values, te1=None, te2=None): # For convienience if both are missing we should give one error about them if te1 is None and te2 is None: - raise RuntimeError('EchoTime1 and EchoTime2 metadata fields not found. ' - 'Please consult the BIDS specification.') + raise RuntimeError( + "EchoTime1 and EchoTime2 metadata fields not found. " + "Please consult the BIDS specification." + ) if te1 is None: raise RuntimeError( - 'EchoTime1 metadata field not found. Please consult the BIDS specification.') + "EchoTime1 metadata field not found. Please consult the BIDS specification." + ) if te2 is None: raise RuntimeError( - 'EchoTime2 metadata field not found. Please consult the BIDS specification.') + "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(caching='unchanged') # Read as float64 for safety - hdr = im.header.copy() - - # Rescale to [0, 2*pi] - data = (data - data.min()) * (2 * np.pi / (data.max() - data.min())) - - # Round to float32 and clip - data = np.clip(np.float32(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, None, hdr).to_filename(out_file) - return out_file - - -def _subtract_phases(in_phases, in_meta, newpath=None): - # Discard traits with copy(), so that pop() works. - in_meta = (in_meta[0].copy(), in_meta[1].copy()) - echo_times = tuple([m.pop('EchoTime', None) for m in in_meta]) - if not all(echo_times): - raise ValueError( - 'One or more missing EchoTime metadata parameter ' - 'associated to one or more phase map(s).') - - if echo_times[0] > echo_times[1]: - in_phases = (in_phases[1], in_phases[0]) - in_meta = (in_meta[1], in_meta[0]) - echo_times = (echo_times[1], echo_times[0]) - - in_phases_nii = [nb.load(ph) for ph in in_phases] - sub_data = in_phases_nii[1].get_fdata(dtype='float32') - \ - in_phases_nii[0].get_fdata(dtype='float32') - - # wrap negative radians back to [0, 2pi] - sub_data[sub_data < 0] += 2 * np.pi - sub_data = np.clip(sub_data, 0.0, 2 * np.pi) - - new_meta = in_meta[1].copy() - new_meta.update(in_meta[0]) - new_meta['EchoTime1'] = echo_times[0] - new_meta['EchoTime2'] = echo_times[1] - - hdr = in_phases_nii[0].header.copy() - hdr.set_data_dtype(np.float32) - hdr.set_xyzt_units('mm') - nii = nb.Nifti1Image(sub_data, in_phases_nii[0].affine, hdr) - out_phdiff = fname_presuffix(in_phases[0], suffix='_phdiff', - newpath=newpath) - nii.to_filename(out_phdiff) - return out_phdiff, new_meta diff --git a/sdcflows/utils/__init__.py b/sdcflows/utils/__init__.py new file mode 100644 index 0000000000..e69de29bb2 diff --git a/sdcflows/utils/phasemanip.py b/sdcflows/utils/phasemanip.py new file mode 100644 index 0000000000..59892258a1 --- /dev/null +++ b/sdcflows/utils/phasemanip.py @@ -0,0 +1,109 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +"""Utilities to manipulate phase and phase difference maps.""" + + +def au2rads(in_file, newpath=None): + """Convert the input phase map in arbitrary units (a.u.) to rads.""" + import numpy as np + import nibabel as nb + from nipype.utils.filemanip import fname_presuffix + + im = nb.load(in_file) + data = im.get_fdata(caching="unchanged") # Read as float64 for safety + hdr = im.header.copy() + + # Rescale to [0, 2*pi] + data = (data - data.min()) * (2 * np.pi / (data.max() - data.min())) + + # Round to float32 and clip + data = np.clip(np.float32(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, None, hdr).to_filename(out_file) + return out_file + + +def subtract_phases(in_phases, in_meta, newpath=None): + """Calculate the phase-difference map, given two input phase maps.""" + import numpy as np + import nibabel as nb + from nipype.utils.filemanip import fname_presuffix + + echo_times = tuple([m.pop("EchoTime", None) for m in in_meta]) + if not all(echo_times): + raise ValueError( + "One or more missing EchoTime metadata parameter " + "associated to one or more phase map(s)." + ) + + if echo_times[0] > echo_times[1]: + in_phases = (in_phases[1], in_phases[0]) + in_meta = (in_meta[1], in_meta[0]) + echo_times = (echo_times[1], echo_times[0]) + + in_phases_nii = [nb.load(ph) for ph in in_phases] + sub_data = in_phases_nii[1].get_fdata(dtype="float32") - in_phases_nii[0].get_fdata( + dtype="float32" + ) + + # wrap negative radians back to [0, 2pi] + sub_data[sub_data < 0] += 2 * np.pi + sub_data = np.clip(sub_data, 0.0, 2 * np.pi) + + new_meta = in_meta[1].copy() + new_meta.update(in_meta[0]) + new_meta["EchoTime1"] = echo_times[0] + new_meta["EchoTime2"] = echo_times[1] + + hdr = in_phases_nii[0].header.copy() + hdr.set_data_dtype(np.float32) + hdr.set_xyzt_units("mm") + nii = nb.Nifti1Image(sub_data, in_phases_nii[0].affine, hdr) + out_phdiff = fname_presuffix(in_phases[0], suffix="_phdiff", newpath=newpath) + nii.to_filename(out_phdiff) + return out_phdiff, new_meta + + +def phdiff2fmap(in_file, delta_te, newpath=None): + r""" + Convert the input phase-difference map into a fieldmap in Hz. + + Uses eq. (1) of [Hutton2002]_: + + .. math:: + + \Delta B_0 (\text{T}^{-1}) = \frac{\Delta \Theta}{2\pi\gamma \Delta\text{TE}} + + + In this case, we do not take into account the gyromagnetic ratio of the + proton (:math:`\gamma`), since it will be applied inside FUGUE: + + .. math:: + + \Delta B_0 (\text{Hz}) = \frac{\Delta \Theta}{2\pi \Delta\text{TE}} + + References + ---------- + .. [Hutton2002] Hutton et al., Image Distortion Correction in fMRI: A Quantitative + Evaluation, NeuroImage 16(1):217-240, 2002. doi:`10.1006/nimg.2001.1054 + `_. + + + """ + import math + import numpy as np + import nibabel as nb + from nipype.utils.filemanip import fname_presuffix + + # GYROMAG_RATIO_H_PROTON_MHZ = 42.576 + + out_file = fname_presuffix(in_file, suffix="_fmap", newpath=newpath) + image = nb.load(in_file) + data = image.get_fdata(dtype="float32") / (2.0 * math.pi * delta_te) + nii = nb.Nifti1Image(data, image.affine, image.header) + nii.set_data_dtype(np.float32) + nii.to_filename(out_file) + return out_file