diff --git a/.circleci/config.yml b/.circleci/config.yml index fedb9e92932..4a44a61b9c7 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -151,8 +151,8 @@ jobs: steps: - restore_cache: keys: - - regression-v4-{{ .Revision }} - - regression-v4- + - regression-v5-{{ .Revision }} + - regression-v5- - run: name: Get truncated BOLD series command: | @@ -175,7 +175,7 @@ jobs: echo "Pre-computed masks were cached" fi - save_cache: - key: regression-v4-{{ .Revision }}-{{ epoch }} + key: regression-v5-{{ .Revision }}-{{ epoch }} paths: - /tmp/data @@ -284,7 +284,7 @@ jobs: - restore_cache: keys: - - regression-v4-{{ .Revision }} + - regression-v5-{{ .Revision }} - restore_cache: keys: - masks-workdir-v2-{{ .Branch }}-{{epoch}} diff --git a/niworkflows/anat/ants.py b/niworkflows/anat/ants.py index c32b3042504..09ec9330246 100644 --- a/niworkflows/anat/ants.py +++ b/niworkflows/anat/ants.py @@ -16,6 +16,7 @@ from nipype.interfaces.ants import N4BiasFieldCorrection, Atropos, MultiplyImages from ..utils.misc import get_template_specs +from ..utils.connections import pop_file as _pop # niworkflows from ..interfaces.ants import ( @@ -897,12 +898,6 @@ def init_n4_only_wf( return wf -def _pop(in_files): - if isinstance(in_files, (list, tuple)): - return in_files[0] - return in_files - - def _select_labels(in_segm, labels): from os import getcwd import numpy as np diff --git a/niworkflows/func/tests/test_util.py b/niworkflows/func/tests/test_util.py index ef59d3f517e..2151fb976a4 100755 --- a/niworkflows/func/tests/test_util.py +++ b/niworkflows/func/tests/test_util.py @@ -1,4 +1,4 @@ -""" Testing module for fmriprep.workflows.bold.util """ +"""Testing module for fmriprep.workflows.bold.util.""" import pytest import os from pathlib import Path @@ -8,10 +8,41 @@ from nipype.utils.filemanip import fname_presuffix, copyfile from nilearn.image import load_img +from ...utils.connections import listify from niworkflows.interfaces.masks import ROIsPlot from ..util import init_bold_reference_wf +# Multi-echo datasets +bold_datasets = ["""\ +ds000210/sub-06_task-rest_run-01_echo-1_bold.nii.gz +ds000210/sub-06_task-rest_run-01_echo-2_bold.nii.gz +ds000210/sub-06_task-rest_run-01_echo-3_bold.nii.gz\ +""".splitlines(), """\ +ds000216/sub-03_task-rest_echo-1_bold.nii.gz +ds000216/sub-03_task-rest_echo-2_bold.nii.gz +ds000216/sub-03_task-rest_echo-3_bold.nii.gz +ds000216/sub-03_task-rest_echo-4_bold.nii.gz""".splitlines()] + +# Single-echo datasets +bold_datasets += """\ +ds000116/sub-12_task-visualoddballwithbuttonresponsetotargetstimuli_run-02_bold.nii.gz +ds000133/sub-06_ses-post_task-rest_run-01_bold.nii.gz +ds000140/sub-32_task-heatpainwithregulationandratings_run-02_bold.nii.gz +ds000157/sub-23_task-passiveimageviewing_bold.nii.gz +ds000237/sub-03_task-MemorySpan_acq-multiband_run-01_bold.nii.gz +ds000237/sub-06_task-MemorySpan_acq-multiband_run-01_bold.nii.gz +ds001240/sub-26_task-localizerimagination_bold.nii.gz +ds001240/sub-26_task-localizerviewing_bold.nii.gz +ds001240/sub-26_task-molencoding_run-01_bold.nii.gz +ds001240/sub-26_task-molencoding_run-02_bold.nii.gz +ds001240/sub-26_task-molretrieval_run-01_bold.nii.gz +ds001240/sub-26_task-molretrieval_run-02_bold.nii.gz +ds001240/sub-26_task-rest_bold.nii.gz +ds001362/sub-01_task-taskname_run-01_bold.nii.gz""".splitlines() + +bold_datasets = [listify(d) for d in bold_datasets] + def symmetric_overlap(img1, img2): mask1 = load_img(img1).get_fdata() > 0 @@ -32,43 +63,23 @@ def symmetric_overlap(img1, img2): "input_fname,expected_fname", [ ( - os.path.join(os.getenv("FMRIPREP_REGRESSION_SOURCE", ""), base_fname), + [os.path.join(os.getenv("FMRIPREP_REGRESSION_SOURCE", ""), bf) + for bf in base_fname], fname_presuffix( - base_fname, + base_fname[0].replace("_echo-1", ""), suffix="_mask", use_ext=True, newpath=os.path.join( os.getenv("FMRIPREP_REGRESSION_TARGETS", ""), - os.path.dirname(base_fname), + os.path.dirname(base_fname[0]), ), ), ) - for base_fname in """\ -ds000116/sub-12_task-visualoddballwithbuttonresponsetotargetstimuli_run-02_bold.nii.gz -ds000133/sub-06_ses-post_task-rest_run-01_bold.nii.gz -ds000140/sub-32_task-heatpainwithregulationandratings_run-02_bold.nii.gz -ds000157/sub-23_task-passiveimageviewing_bold.nii.gz -ds000210/sub-06_task-rest_run-01_echo-1_bold.nii.gz -ds000210/sub-06_task-rest_run-01_echo-2_bold.nii.gz -ds000210/sub-06_task-rest_run-01_echo-3_bold.nii.gz -ds000216/sub-03_task-rest_echo-1_bold.nii.gz -ds000216/sub-03_task-rest_echo-2_bold.nii.gz -ds000216/sub-03_task-rest_echo-3_bold.nii.gz -ds000216/sub-03_task-rest_echo-4_bold.nii.gz -ds000237/sub-03_task-MemorySpan_acq-multiband_run-01_bold.nii.gz -ds000237/sub-06_task-MemorySpan_acq-multiband_run-01_bold.nii.gz -ds001240/sub-26_task-localizerimagination_bold.nii.gz -ds001240/sub-26_task-localizerviewing_bold.nii.gz -ds001240/sub-26_task-molencoding_run-01_bold.nii.gz -ds001240/sub-26_task-molencoding_run-02_bold.nii.gz -ds001240/sub-26_task-molretrieval_run-01_bold.nii.gz -ds001240/sub-26_task-molretrieval_run-02_bold.nii.gz -ds001240/sub-26_task-rest_bold.nii.gz -ds001362/sub-01_task-taskname_run-01_bold.nii.gz""".splitlines() + for base_fname in bold_datasets ], ) def test_masking(input_fname, expected_fname): - basename = Path(input_fname).name + basename = Path(input_fname[0]).name dsname = Path(expected_fname).parent.name # Reconstruct base_fname from above @@ -76,8 +87,10 @@ def test_masking(input_fname, expected_fname): newpath = reports_dir / dsname name = basename.rstrip("_bold.nii.gz").replace("-", "_") - bold_reference_wf = init_bold_reference_wf(omp_nthreads=1, name=name) - bold_reference_wf.inputs.inputnode.bold_file = input_fname + bold_reference_wf = init_bold_reference_wf(omp_nthreads=1, name=name, + multiecho=len(input_fname) > 1) + bold_reference_wf.inputs.inputnode.bold_file = input_fname[0] if len(input_fname) == 1 \ + else input_fname base_dir = os.getenv("CACHED_WORK_DIRECTORY") if base_dir: base_dir = Path(base_dir) / dsname @@ -85,7 +98,7 @@ def test_masking(input_fname, expected_fname): bold_reference_wf.base_dir = str(base_dir) out_fname = fname_presuffix( - basename, suffix="_mask.svg", use_ext=False, newpath=str(newpath) + Path(expected_fname).name, suffix=".svg", use_ext=False, newpath=str(newpath) ) newpath.mkdir(parents=True, exist_ok=True) @@ -117,7 +130,7 @@ def test_masking(input_fname, expected_fname): mask_dir.mkdir(parents=True, exist_ok=True) copyfile( combine_masks.result.outputs.out_file, - fname_presuffix(basename, suffix="_mask", use_ext=True, newpath=str(mask_dir)), + str(mask_dir / Path(expected_fname).name), copy=True, ) diff --git a/niworkflows/func/util.py b/niworkflows/func/util.py index c2c4496268b..cd9140425da 100644 --- a/niworkflows/func/util.py +++ b/niworkflows/func/util.py @@ -20,6 +20,7 @@ from ..interfaces.masks import SimpleShowMaskRPT from ..interfaces.registration import EstimateReferenceImage from ..interfaces.utils import CopyXForm +from ..utils.connections import listify from ..utils.misc import pass_dummy_scans as _pass_dummy_scans @@ -29,8 +30,10 @@ def init_bold_reference_wf( omp_nthreads, bold_file=None, + sbref_files=None, brainmask_thresh=0.85, pre_mask=False, + multiecho=False, name="bold_reference_wf", gen_report=False, ): @@ -51,19 +54,26 @@ def init_bold_reference_wf( Parameters ---------- - omp_nthreads : int + omp_nthreads : :obj:`int` Maximum number of threads an individual process may use - bold_file : str + bold_file : :obj:`str` BOLD series NIfTI file + sbref_files : :obj:`list` or :obj:`bool` + Single band (as opposed to multi band) reference NIfTI file. + If ``True`` is passed, the workflow is built to accommodate SBRefs, + but the input is left undefined (i.e., it is left open for connection) brainmask_thresh: :obj:`float` Lower threshold for the probabilistic brainmask to obtain the final binary mask (default: 0.85). - pre_mask : bool + pre_mask : :obj:`bool` Indicates whether the ``pre_mask`` input will be set (and thus, step 1 should be skipped). - name : str + multiecho : :obj:`bool` + If multiecho data was supplied, data from the first echo + will be selected + name : :obj:`str` Name of workflow (default: ``bold_reference_wf``) - gen_report : bool + gen_report : :obj:`bool` Whether a mask report node should be appended in the end Inputs @@ -104,10 +114,12 @@ def init_bold_reference_wf( """ workflow = Workflow(name=name) - workflow.__desc__ = """\ + workflow.__desc__ = f"""\ First, a reference volume and its skull-stripped version were generated -using a custom methodology of *fMRIPrep*. +{'from the shortest echo of the BOLD run' * multiecho} using a custom +methodology of *fMRIPrep*. """ + inputnode = pe.Node( niu.IdentityInterface( fields=["bold_file", "bold_mask", "dummy_scans", "sbref_file"] @@ -125,6 +137,7 @@ def init_bold_reference_wf( "ref_image_brain", "bold_mask", "validation_report", + "mask_report", ] ), name="outputnode", @@ -134,10 +147,15 @@ def init_bold_reference_wf( if bold_file is not None: inputnode.inputs.bold_file = bold_file - validate = pe.Node(ValidateImage(), name="validate", mem_gb=DEFAULT_MEMORY_MIN_GB) + val_bold = pe.MapNode( + ValidateImage(), + name="val_bold", + mem_gb=DEFAULT_MEMORY_MIN_GB, + iterfield=["in_file"], + ) gen_ref = pe.Node( - EstimateReferenceImage(), name="gen_ref", mem_gb=1 + EstimateReferenceImage(multiecho=multiecho), name="gen_ref", mem_gb=1 ) # OE: 128x128x128x50 * 64 / 8 ~ 900MB. enhance_and_skullstrip_bold_wf = init_enhance_and_skullstrip_bold_wf( brainmask_thresh=brainmask_thresh, @@ -151,23 +169,23 @@ def init_bold_reference_wf( run_without_submitting=True, mem_gb=DEFAULT_MEMORY_MIN_GB, ) + bold_1st = pe.Node(niu.Select(index=[0]), + name="bold_1st", run_without_submitting=True) + validate_1st = pe.Node(niu.Select(index=[0]), + name="validate_1st", run_without_submitting=True) # fmt: off workflow.connect([ + (inputnode, val_bold, [(("bold_file", listify), "in_file")]), (inputnode, enhance_and_skullstrip_bold_wf, [ ("bold_mask", "inputnode.pre_mask"), ]), - (inputnode, validate, [("bold_file", "in_file")]), - (inputnode, gen_ref, [("sbref_file", "sbref_file")]), (inputnode, calc_dummy_scans, [("dummy_scans", "dummy_scans")]), - (validate, gen_ref, [("out_file", "in_file")]), + (val_bold, gen_ref, [("out_file", "in_file")]), (gen_ref, enhance_and_skullstrip_bold_wf, [ ("ref_image", "inputnode.in_file"), ]), - (validate, outputnode, [ - ("out_file", "bold_file"), - ("out_report", "validation_report"), - ]), + (val_bold, bold_1st, [(("out_file", listify), "inlist")]), (gen_ref, calc_dummy_scans, [("n_volumes_to_discard", "algo_dummy_scans")]), (calc_dummy_scans, outputnode, [("skip_vols_num", "skip_vols")]), (gen_ref, outputnode, [ @@ -179,9 +197,39 @@ def init_bold_reference_wf( ("outputnode.mask_file", "bold_mask"), ("outputnode.skull_stripped_file", "ref_image_brain"), ]), + (val_bold, validate_1st, [(("out_report", listify), "inlist")]), + (bold_1st, outputnode, [("out", "bold_file")]), + (validate_1st, outputnode, [("out", "validation_report")]), ]) # fmt: on + if sbref_files: + nsbrefs = 0 + if sbref_files is not True: + # If not boolean, then it is a list-of or pathlike. + inputnode.inputs.sbref_file = sbref_files + nsbrefs = 1 if isinstance(sbref_files, str) else len(sbref_files) + + val_sbref = pe.MapNode( + ValidateImage(), + name="val_sbref", + mem_gb=DEFAULT_MEMORY_MIN_GB, + iterfield=["in_file"], + ) + # fmt: off + workflow.connect([ + (inputnode, val_sbref, [(("sbref_file", listify), "in_file")]), + (val_sbref, gen_ref, [("out_file", "sbref_file")]), + ]) + # fmt: on + + # Edit the boilerplate as the SBRef will be the reference + workflow.__desc__ = f"""\ +First, a reference volume and its skull-stripped version were generated +by aligning and averaging{' the first echo of' * multiecho} +{nsbrefs or ''} single-band references (SBRefs). +""" + if gen_report: mask_reportlet = pe.Node(SimpleShowMaskRPT(), name="mask_reportlet") # fmt: off diff --git a/niworkflows/interfaces/registration.py b/niworkflows/interfaces/registration.py index 60a072c6757..5a71eee4aec 100644 --- a/niworkflows/interfaces/registration.py +++ b/niworkflows/interfaces/registration.py @@ -16,6 +16,7 @@ BaseInterfaceInputSpec, File, SimpleInterface, + InputMultiObject, ) from nipype.interfaces.mixins import reporting from nipype.interfaces import freesurfer as fs @@ -398,14 +399,39 @@ def _post_run_hook(self, runtime): class _EstimateReferenceImageInputSpec(BaseInterfaceInputSpec): - in_file = File(exists=True, mandatory=True, desc="4D EPI file") - sbref_file = File(exists=True, desc="Single band reference image") + in_file = InputMultiObject( + File(exists=True), + mandatory=True, + desc=( + "4D EPI file. If multiple files " + "are provided, they are assumed " + "to represent multiple echoes " + "from the same run." + ), + ) + sbref_file = InputMultiObject( + File(exists=True), + desc=( + "Single band reference image. " + "If multiple files are provided, " + "they are assumed to represent " + "multiple echoes." + ), + ) mc_method = traits.Enum( "AFNI", "FSL", usedefault=True, desc="Which software to use to perform motion correction", ) + multiecho = traits.Bool( + False, + usedefault=True, + desc=( + "If multiecho data was supplied, data from " + "the first echo will be selected." + ), + ) class _EstimateReferenceImageOutputSpec(TraitedSpec): @@ -419,51 +445,82 @@ class _EstimateReferenceImageOutputSpec(TraitedSpec): class EstimateReferenceImage(SimpleInterface): """ - Given an 4D EPI file estimate an optimal reference image that could be later - used for motion estimation and coregistration purposes. If detected uses - T1 saturated volumes (non-steady state) otherwise a median of + Generate a reference 3D map from BOLD and SBRef EPI images for BOLD datasets. + + Given a 4D BOLD file or one or more 3/4D SBRefs, estimate a reference + image for subsequent motion estimation and coregistration steps. + For the case of BOLD datasets, it estimates a number of T1w saturated volumes + (non-steady state at the beginning of the scan) and calculates the median + across them. + Otherwise (SBRefs or detected zero non-steady state frames), a median of of a subset of motion corrected volumes is used. + If the input reference (BOLD or SBRef) is 3D already, it just returns a + copy of the image with the NIfTI header extensions removed. + + LIMITATION: If one wants to extract the reference from several SBRefs + with several echoes each, the first echo should be selected elsewhere + and run this interface in ``multiecho = False`` mode. """ input_spec = _EstimateReferenceImageInputSpec output_spec = _EstimateReferenceImageOutputSpec def _run_interface(self, runtime): - ref_name = self.inputs.in_file - ref_nii = nb.load(ref_name) - n_volumes_to_discard = _get_vols_to_discard(ref_nii) + is_sbref = isdefined(self.inputs.sbref_file) + ref_input = self.inputs.sbref_file if is_sbref else self.inputs.in_file + + if self.inputs.multiecho: + if len(ref_input) < 2: + input_name = "sbref_file" if is_sbref else "in_file" + raise ValueError("Argument 'multiecho' is True but " + f"'{input_name}' has only one element.") + else: + # Select only the first echo (see LIMITATION above for SBRefs) + ref_input = ref_input[:1] + elif not is_sbref and len(ref_input) > 1: + raise ValueError("Input 'in_file' cannot point to more than one file " + "for single-echo BOLD datasets.") + + # Build the nibabel spatial image we will work with + ref_im = [] + for im_i in ref_input: + nib_i = nb.squeeze_image(nb.load(im_i)) + if nib_i.dataobj.ndim == 3: + ref_im.append(nib_i) + elif nib_i.dataobj.ndim == 4: + ref_im += nb.four_to_three(nib_i) + ref_im = nb.squeeze_image(nb.concat_images(ref_im)) + + # Volumes to discard only makes sense with BOLD inputs. + if not is_sbref: + n_volumes_to_discard = _get_vols_to_discard(ref_im) + out_ref_fname = os.path.join(runtime.cwd, "ref_bold.nii.gz") + else: + n_volumes_to_discard = 0 + out_ref_fname = os.path.join(runtime.cwd, "ref_sbref.nii.gz") + # Set interface outputs self._results["n_volumes_to_discard"] = n_volumes_to_discard - - out_ref_fname = os.path.join(runtime.cwd, "ref_bold.nii.gz") - if isdefined(self.inputs.sbref_file): - out_ref_fname = os.path.join(runtime.cwd, "ref_sbref.nii.gz") - ref_name = self.inputs.sbref_file - ref_nii = nb.squeeze_image(nb.load(ref_name)) - - # If reference is only 1 volume, return it directly - if len(ref_nii.shape) == 3: - ref_nii.header.extensions.clear() - ref_nii.to_filename(out_ref_fname) - self._results["ref_image"] = out_ref_fname - return runtime - else: - # Reset this variable as it no longer applies - # and value for the output is stored in self._results - n_volumes_to_discard = 0 + self._results["ref_image"] = out_ref_fname # Slicing may induce inconsistencies with shape-dependent values in extensions. # For now, remove all. If this turns out to be a mistake, we can select extensions # that don't break pipeline stages. - ref_nii.header.extensions.clear() + ref_im.header.extensions.clear() + + # If reference is only 1 volume, return it directly + if ref_im.dataobj.ndim == 3: + ref_im.to_filename(out_ref_fname) + return runtime if n_volumes_to_discard == 0: - if ref_nii.shape[-1] > 40: - ref_name = os.path.join(runtime.cwd, "slice.nii.gz") - nb.Nifti1Image( - ref_nii.dataobj[:, :, :, 20:40], ref_nii.affine, ref_nii.header - ).to_filename(ref_name) + if ref_im.shape[-1] > 40: + ref_im = nb.Nifti1Image( + ref_im.dataobj[:, :, :, 20:40], ref_im.affine, ref_im.header + ) + ref_name = os.path.join(runtime.cwd, "slice.nii.gz") + ref_im.to_filename(ref_name) if self.inputs.mc_method == "AFNI": res = afni.Volreg( in_file=ref_name, @@ -480,15 +537,12 @@ def _run_interface(self, runtime): median_image_data = np.median(mc_slice_nii.get_fdata(), axis=3) else: median_image_data = np.median( - ref_nii.dataobj[:, :, :, :n_volumes_to_discard], axis=3 + ref_im.dataobj[:, :, :, :n_volumes_to_discard], axis=3 ) - nb.Nifti1Image(median_image_data, ref_nii.affine, ref_nii.header).to_filename( + nb.Nifti1Image(median_image_data, ref_im.affine, ref_im.header).to_filename( out_ref_fname ) - - self._results["ref_image"] = out_ref_fname - return runtime diff --git a/niworkflows/utils/connections.py b/niworkflows/utils/connections.py new file mode 100644 index 00000000000..7862cf4fd74 --- /dev/null +++ b/niworkflows/utils/connections.py @@ -0,0 +1,58 @@ +# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*- +# vi: set ft=python sts=4 ts=4 sw=4 et: +""" +Utilities for the creation of nipype workflows. + +Because these functions are meant to be inlined in nipype's ``connect`` invocations, +all the imports MUST be done in each function's context. + +""" + +__all__ = [ + "listify", + "pop_file", +] + + +def pop_file(in_files): + """ + Select the first file from a list of filenames. + + Used to grab the first echo's file when processing + multi-echo data through workflows that only accept + a single file. + + Examples + -------- + >>> pop_file('some/file.nii.gz') + 'some/file.nii.gz' + >>> pop_file(['some/file1.nii.gz', 'some/file2.nii.gz']) + 'some/file1.nii.gz' + + """ + if isinstance(in_files, (list, tuple)): + return in_files[0] + return in_files + + +def listify(value): + """ + Convert to a list (inspired by bids.utils.listify). + + Examples + -------- + >>> listify('some/file.nii.gz') + ['some/file.nii.gz'] + >>> listify((0.1, 0.2)) + [0.1, 0.2] + >>> listify(None) is None + True + + """ + from pathlib import Path + from nipype.interfaces.base import isdefined + if not isdefined(value) or isinstance(value, type(None)): + return value + if isinstance(value, (str, bytes, Path)): + return [str(value)] + return list(value)