ENH: Resample BOLD data to any surface template space using the Connectome Workbench (#3461)

## Changes proposed in this pull request

Currently in fMRIPrep the Connectome Workbench is used to resample BOLD
data to fsLR spaces, as part of the CIFTI workflow.
This PR allows resampling BOLD data to any surface template space with a
similar workflow, as long as (a) the "space-fsLR" spheres exist for the
template (e.g., "tpl-onavg_space-fsLR_hemi-R_den-41k_sphere.surf.gii")
or (b) the template is "fsLR".

In the long run, we might want to generate alternative generic
workflows, e.g., (a) using the [SurfaceTransform
class](nipy/nitransforms#203) of nitransforms by
@Shotgunosine and me, and (b) using FreeSurfer's resampling workflow.

This PR depends on [a recent PR of
smriprep](nipreps/smriprep#473) and updating
TemplateFlow's files on AWS S3.
Currently [the wrong subfolder on S3 was
updated](https://github.com/templateflow/tpl-onavg/actions/runs/14868227699/job/41750227920#step:3:182)
for some unknown reason.

This generic resampling workflow might serve as the first step for the
[generic CIFTI
workflow](#3330).

## Documentation that should be reviewed

I've updated the workflows' docstrings accordingly, but feel free to
make suggestions/edits.

@effigies @oesteban

Author: effigies

fmriprep/data/boilerplate.bib

@@ -365,3 +365,17 @@ @article{patriat_improved_2017
     keywords = {Motion, Correction, Methods, Rs-fMRI},
     pages = {74--82},
 }
+
+@article{onavg,
+	author = {Feilong, Ma and Jiahui, Guo and Gobbini, Maria Ida and Haxby, James V.},
+	title = {A cortical surface template for human neuroscience},
+	url = {https://www.nature.com/articles/s41592-024-02346-y},
+	journal = {Nature Methods},
+	issn = {1548-7105},
+	number = {9},
+	volume = {21},
+	year = {2024},
+	month = sep,
+	pages = {1736--1742},
+	doi = {10.1038/s41592-024-02346-y},
+}

fmriprep/workflows/base.py

@@ -532,9 +532,10 @@ def init_single_subject_wf(
                 grayord_density=config.workflow.cifti_output,
                 omp_nthreads=omp_nthreads,
             )
+            fslr_density = '32k' if config.workflow.cifti_output == '91k' else '59k'
             resample_surfaces_wf = init_resample_surfaces_wf(
                 surfaces=['white', 'pial', 'midthickness'],
-                grayord_density=config.workflow.cifti_output,
+                density=fslr_density,
             )
             ds_grayord_metrics_wf = init_ds_grayord_metrics_wf(
                 bids_root=bids_root,
@@ -547,7 +548,7 @@ def init_single_subject_wf(
                 surfaces=['white', 'pial', 'midthickness'],
                 entities={
                     'space': 'fsLR',
-                    'density': '32k' if config.workflow.cifti_output == '91k' else '59k',
+                    'density': fslr_density,
                 },
                 name='ds_fsLR_surfaces_wf',
             )

fmriprep/workflows/bold/base.py

@@ -279,6 +279,7 @@ def init_bold_wf(
     spaces = config.workflow.spaces
     nonstd_spaces = set(spaces.get_nonstandard())
     freesurfer_spaces = spaces.get_fs_spaces()
+    surf_std = [x for x in spaces.get_standard(dim=(2,)) if x.space != 'fsaverage']
 
     #
     # Resampling outputs workflow:
@@ -503,6 +504,39 @@ def init_bold_wf(
             (merge_bold_sources, ds_bold_std_wf, [('out', 'inputnode.source_files')]),
         ])  # fmt:skip
 
+    # Goodvoxels mask might be needed in any surface resampling
+    if config.workflow.project_goodvoxels and (config.workflow.cifti_output or surf_std):
+        from .resampling import init_goodvoxels_bold_mask_wf
+
+        goodvoxels_bold_mask_wf = init_goodvoxels_bold_mask_wf(mem_gb['resampled'])
+        ds_goodvoxels_mask = pe.Node(
+            DerivativesDataSink(
+                base_directory=fmriprep_dir,
+                dismiss_entities=dismiss_echo(),
+                compress=True,
+                space='T1w',
+                desc='goodvoxels',
+                suffix='mask',
+            ),
+            name='ds_goodvoxels_mask',
+            run_without_submitting=True,
+        )
+        ds_goodvoxels_mask.inputs.source_file = bold_file
+
+        workflow.__postdesc__ += """\
+A "goodvoxels" mask was applied during volume-to-surface sampling, excluding
+voxels whose time-series have a locally high coefficient of variation.
+"""
+        workflow.connect([
+            (inputnode, goodvoxels_bold_mask_wf, [('anat_ribbon', 'inputnode.anat_ribbon')]),
+            (bold_anat_wf, goodvoxels_bold_mask_wf, [
+                ('outputnode.bold_file', 'inputnode.bold_file'),
+            ]),
+            (goodvoxels_bold_mask_wf, ds_goodvoxels_mask, [
+                    ('outputnode.goodvoxels_mask', 'in_file'),
+                ]),
+        ])  # fmt:skip
+
     if config.workflow.run_reconall and freesurfer_spaces:
         workflow.__postdesc__ += """\
 Non-gridded (surface) resamplings were performed using `mri_vol2surf`
@@ -550,7 +584,6 @@ def init_bold_wf(
         from .resampling import (
             init_bold_fsLR_resampling_wf,
             init_bold_grayords_wf,
-            init_goodvoxels_bold_mask_wf,
         )
 
         bold_MNI6_wf = init_bold_volumetric_resample_wf(
@@ -569,42 +602,12 @@ def init_bold_wf(
         )
 
         if config.workflow.project_goodvoxels:
-            goodvoxels_bold_mask_wf = init_goodvoxels_bold_mask_wf(mem_gb['resampled'])
-
-            workflow.connect([
-                (inputnode, goodvoxels_bold_mask_wf, [('anat_ribbon', 'inputnode.anat_ribbon')]),
-                (bold_anat_wf, goodvoxels_bold_mask_wf, [
-                    ('outputnode.bold_file', 'inputnode.bold_file'),
-                ]),
-            ])  # fmt:skip
-
-            ds_goodvoxels_mask = pe.Node(
-                DerivativesDataSink(
-                    base_directory=fmriprep_dir,
-                    dismiss_entities=dismiss_echo(),
-                    compress=True,
-                    space='T1w',
-                    desc='goodvoxels',
-                    suffix='mask',
-                ),
-                name='ds_goodvoxels_mask',
-                run_without_submitting=True,
-            )
-            ds_goodvoxels_mask.inputs.source_file = bold_file
             workflow.connect([
-                (goodvoxels_bold_mask_wf, ds_goodvoxels_mask, [
-                    ('outputnode.goodvoxels_mask', 'in_file'),
-                ]),
                 (goodvoxels_bold_mask_wf, bold_fsLR_resampling_wf, [
                     ('outputnode.goodvoxels_mask', 'inputnode.volume_roi'),
                 ]),
             ])  # fmt:skip
 
-            bold_fsLR_resampling_wf.__desc__ += """\
-A "goodvoxels" mask was applied during volume-to-surface sampling in fsLR space,
-excluding voxels whose time-series have a locally high coefficient of variation.
-"""
-
         bold_grayords_wf = init_bold_grayords_wf(
             grayord_density=config.workflow.cifti_output,
             mem_gb=1,
@@ -670,6 +673,105 @@ def init_bold_wf(
             ]),
         ])  # fmt:skip
 
+        if surf_std:
+            from smriprep.workflows.surfaces import init_resample_surfaces_wf
+
+            from .resampling import (
+                init_wb_surf_surf_wf,
+                init_wb_vol_surf_wf,
+            )
+
+            workflow.__postdesc__ += (
+                'Non-gridded (surface) resamplings were performed using the Connectome Workbench.'
+            )
+            config.loggers.workflow.debug('Creating BOLD surface workbench resampling workflow.')
+
+            wb_vol_surf_wf = init_wb_vol_surf_wf(
+                omp_nthreads=omp_nthreads,
+                mem_gb=mem_gb['resampled'],
+                dilate=True,
+            )
+            workflow.connect([
+                (inputnode, wb_vol_surf_wf,[
+                    ('white', 'inputnode.white'),
+                    ('pial', 'inputnode.pial'),
+                    ('midthickness', 'inputnode.midthickness'),
+                ]),
+                (bold_anat_wf, wb_vol_surf_wf, [
+                    ('outputnode.bold_file', 'inputnode.bold_file'),
+                ]),
+            ])  # fmt:skip
+
+            if config.workflow.project_goodvoxels:
+                workflow.connect([
+                    (goodvoxels_bold_mask_wf, wb_vol_surf_wf, [
+                        ('outputnode.goodvoxels_mask', 'inputnode.volume_roi'),
+                    ]),
+                ])  # fmt:skip
+
+            for ref_ in surf_std:
+                template = ref_.space
+                density = ref_.spec.get('density') or ref_.spec.get('den') or None
+                if density is None:
+                    config.loggers.warning(f'Cannot resample {ref_} without density specified.')
+                    continue
+
+                resample_surfaces_wf = init_resample_surfaces_wf(
+                    name=f'resample_surfaces_wf_{template}_{density}',
+                    surfaces=['midthickness'],
+                    template=template,
+                    density=density,
+                )
+
+                wb_surf_surf_wf = init_wb_surf_surf_wf(
+                    template=template,
+                    density=density,
+                    omp_nthreads=omp_nthreads,
+                    mem_gb=mem_gb['resampled'],
+                )
+
+                ds_bold_surf_wb = pe.Node(
+                    DerivativesDataSink(
+                        base_directory=fmriprep_dir,
+                        hemi=['L', 'R'],
+                        dismiss_entities=dismiss_echo(),
+                        space=template,
+                        density=density,
+                        suffix='bold',
+                        TaskName=all_metadata[0].get('TaskName'),
+                        extension='.func.gii',
+                        **prepare_timing_parameters(all_metadata[0]),
+                    ),
+                    iterfield=('in_file', 'hemi'),
+                    name=f'ds_bold_surf_wb_{template}_{density}',
+                    run_without_submitting=True,
+                )
+                ds_bold_surf_wb.inputs.source_file = bold_file
+
+                workflow.connect([
+                    (inputnode, resample_surfaces_wf, [
+                        ('midthickness', 'inputnode.midthickness'),
+                        ('sphere_reg_fsLR', 'inputnode.sphere_reg_fsLR'),
+                    ]),
+                    (wb_vol_surf_wf, wb_surf_surf_wf, [
+                        ('outputnode.bold_fsnative', 'inputnode.bold_fsnative'),
+                    ]),
+                    (inputnode, wb_surf_surf_wf, [
+                        ('midthickness', 'inputnode.midthickness'),
+                        ('sphere_reg_fsLR', 'inputnode.sphere_reg_fsLR'),
+                    ]),
+                    (resample_surfaces_wf, wb_surf_surf_wf, [
+                        (
+                            f'outputnode.midthickness_{template}',
+                            'inputnode.midthickness_resampled'
+                        ),
+                    ]),
+                    (wb_surf_surf_wf, ds_bold_surf_wb, [
+                        ('outputnode.bold_resampled', 'in_file'),
+                        # TODO: json metadata?
+                    ]),
+                ])  # fmt:skip
+
     bold_confounds_wf = init_bold_confs_wf(
         mem_gb=mem_gb['largemem'],
         metadata=all_metadata[0],