ENH: Restore CIFTI-2 generation (#3129)

Builds on #3126.

This one splits out the selection of the MNI template for the CIFTI-2
generation from the rest of the resampling to volumetric templates. If
someone requests MNI152NLin6Asym and CIFTI, then we will resample twice
in the working directory, but given that resampling is significantly
faster now and does not scale files per-time point, I'm okay with that.
I think it's a pretty rare case, and I suspect we'll want to resample a
smaller ROI than the full brain at some point, so a separate pathway for
that seems okay.

Author: effigies

fmriprep/interfaces/gifti.py

@@ -24,7 +24,7 @@ class CreateROIOutputSpec(TraitedSpec):
 
 
 class CreateROI(SimpleInterface):
-    """Prepare GIFTI shape file for use in"""
+    """Prepare GIFTI thickness file for use as a cortical ROI"""
 
     input_spec = CreateROIInputSpec
     output_spec = CreateROIOutputSpec
@@ -46,17 +46,21 @@ def _run_interface(self, runtime):
         # wb_command -metric-math "abs(var * -1) > 0"
         roi = np.abs(darray.data) > 0
 
+        # Divergence: Set datatype to uint8, since the values are boolean
+        # wb_command sets datatype to float32
         darray = nb.gifti.GiftiDataArray(
             roi,
             intent=darray.intent,
-            datatype=darray.datatype,
+            datatype='uint8',
             encoding=darray.encoding,
             endian=darray.endian,
             coordsys=darray.coordsys,
             ordering=darray.ind_ord,
             meta=meta,
         )
 
+        img.darrays[0] = darray
+
         out_filename = os.path.join(runtime.cwd, f"{subject}.{hemi}.roi.native.shape.gii")
         img.to_filename(out_filename)
         self._results["roi_file"] = out_filename

fmriprep/workflows/base.py

@@ -147,6 +147,7 @@ def init_single_subject_wf(subject_id: str):
     from niworkflows.engine.workflows import LiterateWorkflow as Workflow
     from niworkflows.interfaces.bids import BIDSDataGrabber, BIDSInfo
     from niworkflows.interfaces.nilearn import NILEARN_VERSION
+    from niworkflows.interfaces.utility import KeySelect
     from niworkflows.utils.bids import collect_data
     from niworkflows.utils.misc import fix_multi_T1w_source_name
     from niworkflows.utils.spaces import Reference
@@ -332,8 +333,9 @@ def init_single_subject_wf(subject_id: str):
     ])  # fmt:skip
 
     # Set up the template iterator once, if used
+    template_iterator_wf = None
     if config.workflow.level == "full":
-        if spaces.get_spaces(nonstandard=False, dim=(3,)):
+        if spaces.cached.get_spaces(nonstandard=False, dim=(3,)):
             template_iterator_wf = init_template_iterator_wf(spaces=spaces)
             workflow.connect([
                 (anat_fit_wf, template_iterator_wf, [
@@ -342,6 +344,34 @@ def init_single_subject_wf(subject_id: str):
                 ]),
             ])  # fmt:skip
 
+        # Thread MNI152NLin6Asym standard outputs to CIFTI subworkflow, skipping
+        # the iterator, which targets only output spaces.
+        # This can lead to duplication in the working directory if people actually
+        # want MNI152NLin6Asym outputs, but we'll live with it.
+        if config.workflow.cifti_output:
+            from smriprep.interfaces.templateflow import TemplateFlowSelect
+
+            ref = Reference(
+                "MNI152NLin6Asym",
+                {"res": 2 if config.workflow.cifti_output == "91k" else 1},
+            )
+
+            select_MNI6_xfm = pe.Node(
+                KeySelect(fields=["anat2std_xfm"], key=ref.fullname),
+                name="select_MNI6",
+                run_without_submitting=True,
+            )
+            select_MNI6_tpl = pe.Node(
+                TemplateFlowSelect(template=ref.fullname, resolution=ref.spec['res']),
+                name="select_MNI6_tpl",
+            )
+            workflow.connect([
+                (anat_fit_wf, select_MNI6_xfm, [
+                    ("outputnode.anat2std_xfm", "anat2std_xfm"),
+                    ("outputnode.template", "keys"),
+                ]),
+            ])  # fmt:skip
+
     if config.workflow.anat_only:
         return clean_datasinks(workflow)
 
@@ -362,7 +392,6 @@ def init_single_subject_wf(subject_id: str):
             f"{[e.method for e in fmap_estimators]}."
         )
 
-        from niworkflows.interfaces.utility import KeySelect
         from sdcflows import fieldmaps as fm
         from sdcflows.workflows.base import init_fmap_preproc_wf
 
@@ -507,6 +536,12 @@ def init_single_subject_wf(subject_id: str):
                 ('outputnode.subjects_dir', 'inputnode.subjects_dir'),
                 ('outputnode.subject_id', 'inputnode.subject_id'),
                 ('outputnode.fsnative2t1w_xfm', 'inputnode.fsnative2t1w_xfm'),
+                ('outputnode.white', 'inputnode.white'),
+                ('outputnode.pial', 'inputnode.pial'),
+                ('outputnode.midthickness', 'inputnode.midthickness'),
+                ('outputnode.thickness', 'inputnode.thickness'),
+                ('outputnode.sphere_reg_fsLR', 'inputnode.sphere_reg_fsLR'),
+                ('outputnode.anat_ribbon', 'inputnode.anat_ribbon'),
             ]),
         ])  # fmt:skip
         if fieldmap_id:
@@ -522,16 +557,27 @@ def init_single_subject_wf(subject_id: str):
             ])  # fmt:skip
 
         if config.workflow.level == "full":
-            workflow.connect([
-                (template_iterator_wf, bold_wf, [
-                    ("outputnode.anat2std_xfm", "inputnode.anat2std_xfm"),
-                    ("outputnode.space", "inputnode.std_space"),
-                    ("outputnode.resolution", "inputnode.std_resolution"),
-                    ("outputnode.cohort", "inputnode.std_cohort"),
-                    ("outputnode.std_t1w", "inputnode.std_t1w"),
-                    ("outputnode.std_mask", "inputnode.std_mask"),
-                ]),
-            ])  # fmt:skip
+            if template_iterator_wf is not None:
+                workflow.connect([
+                    (template_iterator_wf, bold_wf, [
+                        ("outputnode.anat2std_xfm", "inputnode.anat2std_xfm"),
+                        ("outputnode.space", "inputnode.std_space"),
+                        ("outputnode.resolution", "inputnode.std_resolution"),
+                        ("outputnode.cohort", "inputnode.std_cohort"),
+                        ("outputnode.std_t1w", "inputnode.std_t1w"),
+                        ("outputnode.std_mask", "inputnode.std_mask"),
+                    ]),
+                ])  # fmt:skip
+
+            # Thread MNI152NLin6Asym standard outputs to CIFTI subworkflow, skipping
+            # the iterator, which targets only output spaces.
+            # This can lead to duplication in the working directory if people actually
+            # want MNI152NLin6Asym outputs, but we'll live with it.
+            if config.workflow.cifti_output:
+                workflow.connect([
+                    (select_MNI6_xfm, bold_wf, [("anat2std_xfm", "inputnode.anat2mni6_xfm")]),
+                    (select_MNI6_tpl, bold_wf, [("brain_mask", "inputnode.mni6_mask")]),
+                ])  # fmt:skip
 
     return clean_datasinks(workflow)
 

fmriprep/workflows/bold/base.py

@@ -54,6 +54,7 @@
     init_ds_registration_wf,
     init_ds_volumes_wf,
     init_func_derivatives_wf,
+    prepare_timing_parameters,
 )
 from .registration import init_bold_reg_wf, init_bold_t1_trans_wf
 from .resampling import (
@@ -189,9 +190,16 @@ def init_bold_wf(
                 "t1w_mask",
                 "t1w_dseg",
                 "t1w_tpms",
+                # FreeSurfer outputs
                 "subjects_dir",
                 "subject_id",
                 "fsnative2t1w_xfm",
+                "white",
+                "midthickness",
+                "pial",
+                "thickness",
+                "sphere_reg_fsLR",
+                "anat_ribbon",
                 # Fieldmap registration
                 "fmap",
                 "fmap_ref",
@@ -206,6 +214,9 @@ def init_bold_wf(
                 "std_cohort",
                 "std_t1w",
                 "std_mask",
+                # MNI152NLin6Asym warp, for CIFTI use
+                "anat2mni6_xfm",
+                "mni6_mask",
             ],
         ),
         name="inputnode",
@@ -389,7 +400,7 @@ def init_bold_wf(
             (bold_anat_wf, ds_bold_t1_wf, [('outputnode.bold_file', 'inputnode.bold')]),
         ])  # fmt:skip
 
-    if spaces.get_spaces(nonstandard=False, dim=(3,)):
+    if spaces.cached.get_spaces(nonstandard=False, dim=(3,)):
         # Missing:
         #  * Clipping BOLD after resampling
         #  * Resampling parcellations
@@ -462,6 +473,87 @@ def init_bold_wf(
             (bold_anat_wf, bold_surf_wf, [("outputnode.bold_file", "inputnode.bold_t1w")]),
         ])  # fmt:skip
 
+    if config.workflow.cifti_output:
+        from .resampling import init_bold_fsLR_resampling_wf, init_bold_grayords_wf
+
+        bold_MNI6_wf = init_bold_volumetric_resample_wf(
+            metadata=all_metadata[0],
+            fieldmap_id=fieldmap_id if not multiecho else None,
+            omp_nthreads=omp_nthreads,
+            name='bold_MNI6_wf',
+        )
+
+        bold_fsLR_resampling_wf = init_bold_fsLR_resampling_wf(
+            estimate_goodvoxels=config.workflow.project_goodvoxels,
+            grayord_density=config.workflow.cifti_output,
+            omp_nthreads=omp_nthreads,
+            mem_gb=mem_gb["resampled"],
+        )
+
+        bold_grayords_wf = init_bold_grayords_wf(
+            grayord_density=config.workflow.cifti_output,
+            mem_gb=mem_gb["resampled"],
+            repetition_time=all_metadata[0]["RepetitionTime"],
+        )
+
+        ds_bold_cifti = pe.Node(
+            DerivativesDataSink(
+                base_directory=fmriprep_dir,
+                space='fsLR',
+                density=config.workflow.cifti_output,
+                suffix='bold',
+                compress=False,
+                TaskName=all_metadata[0].get('TaskName'),
+                **prepare_timing_parameters(all_metadata[0]),
+            ),
+            name='ds_bold_cifti',
+            run_without_submitting=True,
+        )
+        ds_bold_cifti.inputs.source_file = bold_file
+
+        workflow.connect([
+            # Resample BOLD to MNI152NLin6Asym, may duplicate bold_std_wf above
+            (inputnode, bold_MNI6_wf, [
+                ("mni6_mask", "inputnode.target_ref_file"),
+                ("mni6_mask", "inputnode.target_mask"),
+                ("anat2mni6_xfm", "inputnode.anat2std_xfm"),
+                ("fmap_ref", "inputnode.fmap_ref"),
+                ("fmap_coeff", "inputnode.fmap_coeff"),
+                ("fmap_id", "inputnode.fmap_id"),
+            ]),
+            (bold_fit_wf, bold_MNI6_wf, [
+                ("outputnode.coreg_boldref", "inputnode.bold_ref_file"),
+                ("outputnode.boldref2fmap_xfm", "inputnode.boldref2fmap_xfm"),
+                ("outputnode.boldref2anat_xfm", "inputnode.boldref2anat_xfm"),
+            ]),
+            (bold_native_wf, bold_MNI6_wf, [
+                ("outputnode.bold_minimal", "inputnode.bold_file"),
+                ("outputnode.motion_xfm", "inputnode.motion_xfm"),
+            ]),
+            # Resample T1w-space BOLD to fsLR surfaces
+            (inputnode, bold_fsLR_resampling_wf, [
+                ("white", "inputnode.white"),
+                ("pial", "inputnode.pial"),
+                ("midthickness", "inputnode.midthickness"),
+                ("thickness", "inputnode.thickness"),
+                ("sphere_reg_fsLR", "inputnode.sphere_reg_fsLR"),
+                ("anat_ribbon", "inputnode.anat_ribbon"),
+            ]),
+            (bold_anat_wf, bold_fsLR_resampling_wf, [
+                ("outputnode.bold_file", "inputnode.bold_file"),
+            ]),
+            (bold_MNI6_wf, bold_grayords_wf, [
+                ("outputnode.bold_file", "inputnode.bold_std"),
+            ]),
+            (bold_fsLR_resampling_wf, bold_grayords_wf, [
+                ("outputnode.bold_fsLR", "inputnode.bold_fsLR"),
+            ]),
+            (bold_grayords_wf, ds_bold_cifti, [
+                ('outputnode.cifti_bold', 'in_file'),
+                (('outputnode.cifti_metadata', _read_json), 'meta_dict'),
+            ]),
+        ])  # fmt:skip
+
     bold_confounds_wf = init_bold_confs_wf(
         mem_gb=mem_gb["largemem"],
         metadata=all_metadata[0],
@@ -657,7 +749,6 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False):
     spaces = config.workflow.spaces
     fmriprep_dir = str(config.execution.fmriprep_dir)
     freesurfer_spaces = spaces.get_fs_spaces()
-    project_goodvoxels = config.workflow.project_goodvoxels and config.workflow.cifti_output
 
     ref_file = bold_file
     wf_name = _get_wf_name(ref_file, "func_preproc")
@@ -745,52 +836,6 @@ def init_func_preproc_wf(bold_file, has_fieldmap=False):
     # SURFACES ##################################################################################
 
     # CIFTI output
-    if config.workflow.cifti_output:
-        from .resampling import init_bold_fsLR_resampling_wf, init_bold_grayords_wf
-
-        bold_fsLR_resampling_wf = init_bold_fsLR_resampling_wf(
-            estimate_goodvoxels=project_goodvoxels,
-            grayord_density=config.workflow.cifti_output,
-            omp_nthreads=omp_nthreads,
-            mem_gb=mem_gb["resampled"],
-        )
-
-        bold_grayords_wf = init_bold_grayords_wf(
-            grayord_density=config.workflow.cifti_output,
-            mem_gb=mem_gb["resampled"],
-            repetition_time=metadata["RepetitionTime"],
-        )
-
-        # fmt:off
-        workflow.connect([
-            (inputnode, bold_fsLR_resampling_wf, [
-                ("surfaces", "inputnode.surfaces"),
-                ("morphometrics", "inputnode.morphometrics"),
-                ("sphere_reg_fsLR", "inputnode.sphere_reg_fsLR"),
-                ("anat_ribbon", "inputnode.anat_ribbon"),
-            ]),
-            (bold_t1_trans_wf, bold_fsLR_resampling_wf, [
-                ("outputnode.bold_t1", "inputnode.bold_file"),
-            ]),
-            (bold_std_trans_wf, bold_grayords_wf, [
-                ("outputnode.bold_std", "inputnode.bold_std"),
-                ("outputnode.spatial_reference", "inputnode.spatial_reference"),
-            ]),
-            (bold_fsLR_resampling_wf, bold_grayords_wf, [
-                ("outputnode.bold_fsLR", "inputnode.bold_fsLR"),
-            ]),
-            (bold_fsLR_resampling_wf, func_derivatives_wf, [
-                ("outputnode.goodvoxels_mask", "inputnode.goodvoxels_mask"),
-            ]),
-            (bold_fsLR_resampling_wf, outputnode, [
-                ("outputnode.weights_text", "weights_text"),
-            ]),
-            (bold_grayords_wf, outputnode, [
-                ("outputnode.cifti_bold", "bold_cifti"),
-                ("outputnode.cifti_metadata", "cifti_metadata"),
-            ]),
-        ])
-        # fmt:on
 
     if spaces.get_spaces(nonstandard=False, dim=(3,)):
         carpetplot_wf = init_carpetplot_wf(
@@ -950,3 +995,10 @@ def get_img_orientation(imgf):
     """Return the image orientation as a string"""
     img = nb.load(imgf)
     return "".join(nb.aff2axcodes(img.affine))
+
+
+def _read_json(in_file):
+    from json import loads
+    from pathlib import Path
+
+    return loads(Path(in_file).read_text())