ENH: Resample BOLD data to any surface template space using the Connectome Workbench

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],