diff --git a/nipype/interfaces/fsl/epi.py b/nipype/interfaces/fsl/epi.py index e7f3ff4318..eeab08371e 100644 --- a/nipype/interfaces/fsl/epi.py +++ b/nipype/interfaces/fsl/epi.py @@ -10,7 +10,7 @@ import nibabel as nb import warnings -from ...utils.filemanip import split_filename +from ...utils.filemanip import split_filename, fname_presuffix from ...utils import NUMPY_MMAP from ..base import traits, TraitedSpec, InputMultiPath, File, isdefined @@ -553,7 +553,7 @@ class EddyInputSpec(FSLCommandInputSpec): exists=True, mandatory=True, argstr="--imain=%s", - desc=("File containing all the images to estimate " "distortions for"), + desc="File containing all the images to estimate distortions for", ) in_mask = File( exists=True, mandatory=True, argstr="--mask=%s", desc="Mask to indicate brain" @@ -562,10 +562,8 @@ class EddyInputSpec(FSLCommandInputSpec): exists=True, mandatory=True, argstr="--index=%s", - desc=( - "File containing indices for all volumes in --imain " - "into --acqp and --topup" - ), + desc="File containing indices for all volumes in --imain " + "into --acqp and --topup", ) in_acqp = File( exists=True, @@ -577,121 +575,246 @@ class EddyInputSpec(FSLCommandInputSpec): exists=True, mandatory=True, argstr="--bvecs=%s", - desc=("File containing the b-vectors for all volumes in " "--imain"), + desc="File containing the b-vectors for all volumes in --imain", ) in_bval = File( exists=True, mandatory=True, argstr="--bvals=%s", - desc=("File containing the b-values for all volumes in " "--imain"), + desc="File containing the b-values for all volumes in --imain", ) out_base = traits.Str( - "eddy_corrected", - argstr="--out=%s", + default_value="eddy_corrected", usedefault=True, - desc=("basename for output (warped) image"), + argstr="--out=%s", + desc="Basename for output image", ) session = File( exists=True, argstr="--session=%s", - desc=("File containing session indices for all volumes in " "--imain"), + desc="File containing session indices for all volumes in --imain", ) in_topup_fieldcoef = File( exists=True, argstr="--topup=%s", requires=["in_topup_movpar"], - desc=("topup file containing the field " "coefficients"), + desc="Topup results file containing the field coefficients", ) in_topup_movpar = File( - exists=True, requires=["in_topup_fieldcoef"], desc="topup movpar.txt file" + exists=True, + requires=["in_topup_fieldcoef"], + desc="Topup results file containing the movement parameters (movpar.txt)", + ) + field = File( + exists=True, + argstr="--field=%s", + desc="Non-topup derived fieldmap scaled in Hz", + ) + field_mat = File( + exists=True, + argstr="--field_mat=%s", + desc="Matrix specifying the relative positions of the fieldmap, " + "--field, and the first volume of the input file, --imain", ) flm = traits.Enum( - "linear", "quadratic", "cubic", argstr="--flm=%s", desc="First level EC model" + "quadratic", + "linear", + "cubic", + usedefault=True, + argstr="--flm=%s", + desc="First level EC model", ) - slm = traits.Enum( - "none", "linear", "quadratic", argstr="--slm=%s", desc="Second level EC model" + "none", + "linear", + "quadratic", + usedefault=True, + argstr="--slm=%s", + desc="Second level EC model", ) - fep = traits.Bool( False, argstr="--fep", desc="Fill empty planes in x- or y-directions" ) - + initrand = traits.Bool( + False, + argstr="--initrand", + desc="Resets rand for when selecting voxels", + min_ver="5.0.10", + ) interp = traits.Enum( "spline", "trilinear", + usedefault=True, argstr="--interp=%s", desc="Interpolation model for estimation step", ) - nvoxhp = traits.Int( - 1000, + default_value=1000, usedefault=True, argstr="--nvoxhp=%s", - desc=("# of voxels used to estimate the " "hyperparameters"), + desc="# of voxels used to estimate the hyperparameters", ) - fudge_factor = traits.Float( - 10.0, + default_value=10.0, usedefault=True, argstr="--ff=%s", - desc=("Fudge factor for hyperparameter " "error variance"), + desc="Fudge factor for hyperparameter error variance", ) - dont_sep_offs_move = traits.Bool( False, argstr="--dont_sep_offs_move", - desc=("Do NOT attempt to separate " "field offset from subject " "movement"), + desc="Do NOT attempt to separate field offset from subject movement", ) - dont_peas = traits.Bool( False, argstr="--dont_peas", - desc="Do NOT perform a post-eddy alignment of " "shells", + desc="Do NOT perform a post-eddy alignment of shells", ) - fwhm = traits.Float( - desc=("FWHM for conditioning filter when estimating " "the parameters"), + desc="FWHM for conditioning filter when estimating the parameters", argstr="--fwhm=%s", ) - niter = traits.Int( 5, usedefault=True, argstr="--niter=%s", desc="Number of iterations" ) - method = traits.Enum( "jac", "lsr", + usedefault=True, argstr="--resamp=%s", - desc=("Final resampling method (jacobian/least " "squares)"), + desc="Final resampling method (jacobian/least squares)", ) + repol = traits.Bool( False, argstr="--repol", desc="Detect and replace outlier slices" ) + outlier_nstd = traits.Int( + argstr="--ol_nstd", + desc="Number of std off to qualify as outlier", + requires=["repol"], + min_ver="5.0.10", + ) + outlier_nvox = traits.Int( + argstr="--ol_nvox", + desc="Min # of voxels in a slice for inclusion in outlier detection", + requires=["repol"], + min_ver="5.0.10", + ) + outlier_type = traits.Enum( + "sw", + "gw", + "both", + argstr="--ol_type", + desc="Type of outliers, slicewise (sw), groupwise (gw) or both (both)", + requires=["repol"], + min_ver="5.0.10", + ) + outlier_pos = traits.Bool( + False, + argstr="--ol_pos", + desc="Consider both positive and negative outliers if set", + requires=["repol"], + min_ver="5.0.10", + ) + outlier_sqr = traits.Bool( + False, + argstr="--ol_sqr", + desc="Consider outliers among sums-of-squared differences if set", + requires=["repol"], + min_ver="5.0.10", + ) + multiband_factor = traits.Int( + argstr="--mb=%s", desc="Multi-band factor", min_ver="5.0.10" + ) + multiband_offset = traits.Enum( + 0, + 1, + -1, + argstr="--mb_offs=%d", + desc="Multi-band offset (-1 if bottom slice removed, 1 if top slice removed", + requires=["multiband_factor"], + min_ver="5.0.10", + ) + + mporder = traits.Int( + argstr="--mporder=%s", + desc="Order of slice-to-vol movement model", + requires=["use_cuda"], + min_ver="5.0.11", + ) + slice2vol_niter = traits.Int( + argstr="--s2v_niter=%d", + desc="Number of iterations for slice-to-vol", + requires=["mporder"], + min_ver="5.0.11", + ) + slice2vol_lambda = traits.Int( + argstr="--s2v_lambda=%d", + desc="Regularisation weight for slice-to-vol movement (reasonable range 1-10)", + requires=["mporder"], + min_ver="5.0.11", + ) + slice2vol_interp = traits.Enum( + "trilinear", + "spline", + argstr="--s2v_interp=%s", + desc="Slice-to-vol interpolation model for estimation step", + requires=["mporder"], + min_ver="5.0.11", + ) + slice_order = traits.File( + exists=True, + argstr="--slspec=%s", + desc="Name of text file completely specifying slice/group acquisition", + requires=["mporder"], + xor=["json"], + min_ver="5.0.11", + ) + json = traits.File( + exists=True, + argstr="--json=%s", + desc="Name of .json text file with information about slice timing", + requires=["mporder"], + xor=["slice_order"], + min_ver="6.0.1", + ) + + estimate_move_by_susceptibility = traits.Bool( + False, + argstr="--estimate_move_by_susceptibility", + desc="Estimate how susceptibility field changes with subject movement", + min_ver="6.0.1", + ) + mbs_niter = traits.Int( + argstr="--mbs_niter=%s", + desc="Number of iterations for MBS estimation", + requires=["estimate_move_by_susceptibility"], + min_ver="6.0.1", + ) + mbs_lambda = traits.Int( + argstr="--mbs_lambda=%s", + desc="Weighting of regularisation for MBS estimation", + requires=["estimate_move_by_susceptibility"], + min_ver="6.0.1", + ) + mbs_ksp = traits.Int( + argstr="--mbs_ksp=%smm", + desc="Knot-spacing for MBS field estimation", + requires=["estimate_move_by_susceptibility"], + min_ver="6.0.1", + ) + num_threads = traits.Int( 1, usedefault=True, nohash=True, desc="Number of openmp threads to use" ) is_shelled = traits.Bool( False, argstr="--data_is_shelled", - desc="Override internal check to ensure that " - "date are acquired on a set of b-value " - "shells", - ) - field = traits.Str( - argstr="--field=%s", - desc="NonTOPUP fieldmap scaled in Hz - filename has " - "to be provided without an extension. TOPUP is " - "strongly recommended", - ) - field_mat = File( - exists=True, - argstr="--field_mat=%s", - desc="Matrix that specifies the relative locations of " - "the field specified by --field and first volume " - "in file --imain", + desc="Override internal check to ensure that date are acquired " + "on a set of b-value shells", ) + use_cuda = traits.Bool(False, desc="Run eddy using cuda gpu") cnr_maps = traits.Bool( False, desc="Output CNR-Maps", argstr="--cnr_maps", min_ver="5.0.10" @@ -707,38 +830,71 @@ class EddyOutputSpec(TraitedSpec): ) out_parameter = File( exists=True, - desc=( - "text file with parameters definining the field and" - "movement for each scan" - ), + desc="Text file with parameters defining the field and movement for each scan", ) out_rotated_bvecs = File( exists=True, desc="File containing rotated b-values for all volumes" ) out_movement_rms = File( - exists=True, desc='Summary of the "total movement" in each volume' + exists=True, desc="Summary of the 'total movement' in each volume" ) out_restricted_movement_rms = File( exists=True, - desc=( - 'Summary of the "total movement" in each volume ' - "disregarding translation in the PE direction" - ), + desc="Summary of the 'total movement' in each volume " + "disregarding translation in the PE direction", ) out_shell_alignment_parameters = File( exists=True, - desc=( - "File containing rigid body movement parameters " - "between the different shells as estimated by a " - "post-hoc mutual information based registration" - ), + desc="Text file containing rigid body movement parameters " + "between the different shells as estimated by a " + "post-hoc mutual information based registration", + ) + out_shell_pe_translation_parameters = File( + exists=True, + desc="Text file containing translation along the PE-direction " + "between the different shells as estimated by a " + "post-hoc mutual information based registration", + ) + out_shell_pe_translation_parameters = File( + exists=True, + desc="Text file containing translation along the PE-direction " + "between the different shells as estimated by a " + "post-hoc mutual information based registration", + ) + out_outlier_map = File( + exists=True, + desc="Matrix where rows represent volumes and columns represent " + 'slices. "0" indicates that scan-slice is not an outlier ' + 'and "1" indicates that it is', + ) + out_outlier_n_stdev_map = File( + exists=True, + desc="Matrix where rows represent volumes and columns represent " + "slices. Values indicate number of standard deviations off the " + "mean difference between observation and prediction is", + ) + out_outlier_n_sqr_stdev_map = File( + exists=True, + desc="Matrix where rows represent volumes and columns represent " + "slices. Values indicate number of standard deivations off the " + "square root of the mean squared difference between observation " + "and prediction is", ) out_outlier_report = File( exists=True, - desc=( - "Text-file with a plain language report on what " - "outlier slices eddy has found" - ), + desc="Text file with a plain language report on what " + "outlier slices eddy has found", + ) + out_outlier_free = File( + exists=True, + desc="4D image file not corrected for susceptibility or eddy-" + "current distortions or subject movement but with outlier " + "slices replaced", + ) + out_movement_over_time = File( + exists=True, + desc="Text file containing translations (mm) and rotations " + "(radians) for each excitation", ) out_cnr_maps = File(exists=True, desc="path/name of file with the cnr_maps") out_residuals = File(exists=True, desc="path/name of file with the residuals") @@ -748,14 +904,16 @@ class Eddy(FSLCommand): """ Interface for FSL eddy, a tool for estimating and correcting eddy currents induced distortions. `User guide - `_ and + `_ and `more info regarding acqp file - `_. + `_. Examples -------- >>> from nipype.interfaces.fsl import Eddy + + Running eddy on a CPU using OpenMP: >>> eddy = Eddy() >>> eddy.inputs.in_file = 'epi.nii' >>> eddy.inputs.in_mask = 'epi_mask.nii' @@ -763,17 +921,36 @@ class Eddy(FSLCommand): >>> eddy.inputs.in_acqp = 'epi_acqp.txt' >>> eddy.inputs.in_bvec = 'bvecs.scheme' >>> eddy.inputs.in_bval = 'bvals.scheme' + >>> eddy.cmdline # doctest: +ELLIPSIS + 'eddy_openmp --flm=quadratic --ff=10.0 \ +--acqp=epi_acqp.txt --bvals=bvals.scheme --bvecs=bvecs.scheme \ +--imain=epi.nii --index=epi_index.txt --mask=epi_mask.nii \ +--interp=spline --resamp=jac --niter=5 --nvoxhp=1000 \ +--out=.../eddy_corrected --slm=none' + + Running eddy on an Nvidia GPU using cuda: >>> eddy.inputs.use_cuda = True >>> eddy.cmdline # doctest: +ELLIPSIS - 'eddy_cuda --ff=10.0 --acqp=epi_acqp.txt --bvals=bvals.scheme \ ---bvecs=bvecs.scheme --imain=epi.nii --index=epi_index.txt \ ---mask=epi_mask.nii --niter=5 --nvoxhp=1000 --out=.../eddy_corrected' - >>> eddy.inputs.use_cuda = False - >>> eddy.cmdline # doctest: +ELLIPSIS - 'eddy_openmp --ff=10.0 --acqp=epi_acqp.txt --bvals=bvals.scheme \ ---bvecs=bvecs.scheme --imain=epi.nii --index=epi_index.txt \ ---mask=epi_mask.nii --niter=5 --nvoxhp=1000 --out=.../eddy_corrected' - >>> res = eddy.run() # doctest: +SKIP + 'eddy_cuda --flm=quadratic --ff=10.0 \ +--acqp=epi_acqp.txt --bvals=bvals.scheme --bvecs=bvecs.scheme \ +--imain=epi.nii --index=epi_index.txt --mask=epi_mask.nii \ +--interp=spline --resamp=jac --niter=5 --nvoxhp=1000 \ +--out=.../eddy_corrected --slm=none' + + Running eddy with slice-to-volume motion correction: + >>> eddy.inputs.mporder = 6 + >>> eddy.inputs.slice2vol_niter = 5 + >>> eddy.inputs.slice2vol_lambda = 1 + >>> eddy.inputs.slice2vol_interp = 'trilinear' + >>> eddy.inputs.slice_order = 'epi_slspec.txt' + >>> eddy.cmdline # doctest: +ELLIPSIS + 'eddy_cuda --flm=quadratic --ff=10.0 \ +--acqp=epi_acqp.txt --bvals=bvals.scheme --bvecs=bvecs.scheme \ +--imain=epi.nii --index=epi_index.txt --mask=epi_mask.nii \ +--interp=spline --resamp=jac --mporder=6 --niter=5 --nvoxhp=1000 \ +--out=.../eddy_corrected --s2v_interp=trilinear --s2v_lambda=1 \ +--s2v_niter=5 --slspec=epi_slspec.txt --slm=none' + >>> res = eddy.run() # doctest: +SKIP """ @@ -826,6 +1003,8 @@ def _run_interface(self, runtime): def _format_arg(self, name, spec, value): if name == "in_topup_fieldcoef": return spec.argstr % value.split("_fieldcoef")[0] + if name == "field": + return spec.argstr % fname_presuffix(value, use_ext=False) if name == "out_base": return spec.argstr % os.path.abspath(value) return super(Eddy, self)._format_arg(name, spec, value) @@ -850,9 +1029,31 @@ def _list_outputs(self): out_shell_alignment_parameters = os.path.abspath( "%s.eddy_post_eddy_shell_alignment_parameters" % self.inputs.out_base ) + out_shell_pe_translation_parameters = os.path.abspath( + "%s.eddy_post_eddy_shell_PE_translation_parameters" % self.inputs.out_base + ) + out_outlier_map = os.path.abspath("%s.eddy_outlier_map" % self.inputs.out_base) + out_outlier_n_stdev_map = os.path.abspath( + "%s.eddy_outlier_n_stdev_map" % self.inputs.out_base + ) + out_outlier_n_sqr_stdev_map = os.path.abspath( + "%s.eddy_outlier_n_sqr_stdev_map" % self.inputs.out_base + ) out_outlier_report = os.path.abspath( "%s.eddy_outlier_report" % self.inputs.out_base ) + if isdefined(self.inputs.repol) and self.inputs.repol: + out_outlier_free = os.path.abspath( + "%s.eddy_outlier_free_data" % self.inputs.out_base + ) + if os.path.exists(out_outlier_free): + outputs["out_outlier_free"] = out_outlier_free + if isdefined(self.inputs.mporder) and self.inputs.mporder > 0: + out_movement_over_time = os.path.abspath( + "%s.eddy_movement_over_time" % self.inputs.out_base + ) + if os.path.exists(out_movement_over_time): + outputs["out_movement_over_time"] = out_movement_over_time if isdefined(self.inputs.cnr_maps) and self.inputs.cnr_maps: out_cnr_maps = os.path.abspath( "%s.eddy_cnr_maps.nii.gz" % self.inputs.out_base @@ -874,6 +1075,16 @@ def _list_outputs(self): outputs["out_restricted_movement_rms"] = out_restricted_movement_rms if os.path.exists(out_shell_alignment_parameters): outputs["out_shell_alignment_parameters"] = out_shell_alignment_parameters + if os.path.exists(out_shell_pe_translation_parameters): + outputs[ + "out_shell_pe_translation_parameters" + ] = out_shell_pe_translation_parameters + if os.path.exists(out_outlier_map): + outputs["out_outlier_map"] = out_outlier_map + if os.path.exists(out_outlier_n_stdev_map): + outputs["out_outlier_n_stdev_map"] = out_outlier_n_stdev_map + if os.path.exists(out_outlier_n_sqr_stdev_map): + outputs["out_outlier_n_sqr_stdev_map"] = out_outlier_n_sqr_stdev_map if os.path.exists(out_outlier_report): outputs["out_outlier_report"] = out_outlier_report diff --git a/nipype/interfaces/fsl/tests/test_auto_Eddy.py b/nipype/interfaces/fsl/tests/test_auto_Eddy.py index cc7eff7a27..4a2d245a23 100644 --- a/nipype/interfaces/fsl/tests/test_auto_Eddy.py +++ b/nipype/interfaces/fsl/tests/test_auto_Eddy.py @@ -9,10 +9,13 @@ def test_Eddy_inputs(): dont_peas=dict(argstr="--dont_peas",), dont_sep_offs_move=dict(argstr="--dont_sep_offs_move",), environ=dict(nohash=True, usedefault=True,), + estimate_move_by_susceptibility=dict( + argstr="--estimate_move_by_susceptibility", min_ver="6.0.1", + ), fep=dict(argstr="--fep",), - field=dict(argstr="--field=%s",), + field=dict(argstr="--field=%s", extensions=None,), field_mat=dict(argstr="--field_mat=%s", extensions=None,), - flm=dict(argstr="--flm=%s",), + flm=dict(argstr="--flm=%s", usedefault=True,), fudge_factor=dict(argstr="--ff=%s", usedefault=True,), fwhm=dict(argstr="--fwhm=%s",), in_acqp=dict(argstr="--acqp=%s", extensions=None, mandatory=True,), @@ -25,18 +28,62 @@ def test_Eddy_inputs(): argstr="--topup=%s", extensions=None, requires=["in_topup_movpar"], ), in_topup_movpar=dict(extensions=None, requires=["in_topup_fieldcoef"],), - interp=dict(argstr="--interp=%s",), + initrand=dict(argstr="--initrand", min_ver="5.0.10",), + interp=dict(argstr="--interp=%s", usedefault=True,), is_shelled=dict(argstr="--data_is_shelled",), - method=dict(argstr="--resamp=%s",), + json=dict( + argstr="--json=%s", + min_ver="6.0.1", + requires=["mporder"], + xor=["slice_order"], + ), + mbs_ksp=dict( + argstr="--mbs_ksp=%smm", + min_ver="6.0.1", + requires=["estimate_move_by_susceptibility"], + ), + mbs_lambda=dict( + argstr="--mbs_lambda=%s", + min_ver="6.0.1", + requires=["estimate_move_by_susceptibility"], + ), + mbs_niter=dict( + argstr="--mbs_niter=%s", + min_ver="6.0.1", + requires=["estimate_move_by_susceptibility"], + ), + method=dict(argstr="--resamp=%s", usedefault=True,), + mporder=dict(argstr="--mporder=%s", min_ver="5.0.11", requires=["use_cuda"],), + multiband_factor=dict(argstr="--mb=%s", min_ver="5.0.10",), + multiband_offset=dict( + argstr="--mb_offs=%d", min_ver="5.0.10", requires=["multiband_factor"], + ), niter=dict(argstr="--niter=%s", usedefault=True,), num_threads=dict(nohash=True, usedefault=True,), nvoxhp=dict(argstr="--nvoxhp=%s", usedefault=True,), out_base=dict(argstr="--out=%s", usedefault=True,), + outlier_nstd=dict(argstr="--ol_nstd", min_ver="5.0.10", requires=["repol"],), + outlier_nvox=dict(argstr="--ol_nvox", min_ver="5.0.10", requires=["repol"],), + outlier_pos=dict(argstr="--ol_pos", min_ver="5.0.10", requires=["repol"],), + outlier_sqr=dict(argstr="--ol_sqr", min_ver="5.0.10", requires=["repol"],), + outlier_type=dict(argstr="--ol_type", min_ver="5.0.10", requires=["repol"],), output_type=dict(), repol=dict(argstr="--repol",), residuals=dict(argstr="--residuals", min_ver="5.0.10",), session=dict(argstr="--session=%s", extensions=None,), - slm=dict(argstr="--slm=%s",), + slice2vol_interp=dict( + argstr="--s2v_interp=%s", min_ver="5.0.11", requires=["mporder"], + ), + slice2vol_lambda=dict( + argstr="--s2v_lambda=%d", min_ver="5.0.11", requires=["mporder"], + ), + slice2vol_niter=dict( + argstr="--s2v_niter=%d", min_ver="5.0.11", requires=["mporder"], + ), + slice_order=dict( + argstr="--slspec=%s", min_ver="5.0.11", requires=["mporder"], xor=["json"], + ), + slm=dict(argstr="--slm=%s", usedefault=True,), use_cuda=dict(), ) inputs = Eddy.input_spec() @@ -50,13 +97,19 @@ def test_Eddy_outputs(): output_map = dict( out_cnr_maps=dict(extensions=None,), out_corrected=dict(extensions=None,), + out_movement_over_time=dict(extensions=None,), out_movement_rms=dict(extensions=None,), + out_outlier_free=dict(extensions=None,), + out_outlier_map=dict(extensions=None,), + out_outlier_n_sqr_stdev_map=dict(extensions=None,), + out_outlier_n_stdev_map=dict(extensions=None,), out_outlier_report=dict(extensions=None,), out_parameter=dict(extensions=None,), out_residuals=dict(extensions=None,), out_restricted_movement_rms=dict(extensions=None,), out_rotated_bvecs=dict(extensions=None,), out_shell_alignment_parameters=dict(extensions=None,), + out_shell_pe_translation_parameters=dict(extensions=None,), ) outputs = Eddy.output_spec() diff --git a/nipype/testing/data/epi_slspec.txt b/nipype/testing/data/epi_slspec.txt new file mode 100644 index 0000000000..e69de29bb2