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

@@ -184,8 +184,9 @@ def init_bold_wf(
         return
 
     config.loggers.workflow.debug(
-        f'Creating bold processing workflow for <{bold_file}> ({mem_gb["filesize"]:.2f} GB / {nvols} TRs). '
-        f'Memory resampled/largemem={mem_gb["resampled"]:.2f}/{mem_gb["largemem"]:.2f} GB.'
+        f'Creating bold processing workflow for <{bold_file}> ({mem_gb["filesize"]:.2f} GB / '
+        f'{nvols} TRs). Memory resampled/largemem={mem_gb["resampled"]:.2f}/'
+        f'{mem_gb["largemem"]:.2f} GB.'
     )
 
     workflow = Workflow(name=_get_wf_name(bold_file, 'bold'))
@@ -494,6 +495,138 @@ def init_bold_wf(
             ]),
         ])  # fmt:skip
 
+    # Goodvoxels mask might be needed in any surface resampling
+    if config.workflow.project_goodvoxels:
+        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
+
+    surf_std = spaces.get_nonstandard(dim=(2,))
+    if surf_std and config.workflow.run_reconall and config.workflow.cifti_output:
+        workflow.__postdesc__ += """\
+Non-gridded (surface) resamplings were performed using the Connectome
+Workbench.
+"""
+        config.loggers.workflow.debug('Creating BOLD surface workbench resampling workflow.')
+        from smriprep.workflows.surfaces import init_resample_surfaces_wf
+
+        from .resampling import (
+            init_wb_surf_surf_wf,
+            init_wb_vol_surf_wf,
+        )
+
+        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:
+            goodvoxels_bold_mask_wf = init_goodvoxels_bold_mask_wf(mem_gb['resampled'])
+
+            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_wb_wf = init_resample_surfaces_wf(
+                name=f'resample_surfaces_wb_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_wb_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_wb_wf, wb_surf_surf_wf, [
+                    ('outputnode.midthickness', 'inputnode.midthickness_resampled'),
+                ]),
+                (wb_surf_surf_wf, ds_bold_surf_wb, [
+                    ('outputnode.bold_resampled', 'in_file'),
+                    # TODO: json metadata?
+                ]),
+            ])  # fmt:skip
+
     if config.workflow.run_reconall and freesurfer_spaces:
         workflow.__postdesc__ += """\
 Non-gridded (surface) resamplings were performed using `mri_vol2surf`
@@ -541,7 +674,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(
@@ -560,42 +692,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,