Merge pull request #408 from tsalo/sbref

ENH: Allow bold_reference_wf to accept multiple EPIs/SBRefs

Author: oesteban

.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}}

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

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,60 +63,42 @@ 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
     reports_dir = Path(os.getenv("FMRIPREP_REGRESSION_REPORTS", ""))
     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
         base_dir.mkdir(parents=True, exist_ok=True)
         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,
     )
 

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