Merge pull request #460 from effigies/enh/fslr_surfs

feat: Output fsLR meshes on subject surfaces

Author: effigies

.circleci/ds005_outputs.txt

@@ -25,6 +25,9 @@ smriprep/sub-01/anat/sub-01_hemi-L_curv.shape.gii
 smriprep/sub-01/anat/sub-01_hemi-L_inflated.surf.gii
 smriprep/sub-01/anat/sub-01_hemi-L_midthickness.surf.gii
 smriprep/sub-01/anat/sub-01_hemi-L_pial.surf.gii
+smriprep/sub-01/anat/sub-01_hemi-L_space-fsLR_den-32k_midthickness.surf.gii
+smriprep/sub-01/anat/sub-01_hemi-L_space-fsLR_den-32k_pial.surf.gii
+smriprep/sub-01/anat/sub-01_hemi-L_space-fsLR_den-32k_white.surf.gii
 smriprep/sub-01/anat/sub-01_hemi-L_space-fsLR_desc-msmsulc_sphere.surf.gii
 smriprep/sub-01/anat/sub-01_hemi-L_space-fsLR_desc-reg_sphere.surf.gii
 smriprep/sub-01/anat/sub-01_hemi-L_space-fsaverage_desc-reg_sphere.surf.gii
@@ -36,6 +39,9 @@ smriprep/sub-01/anat/sub-01_hemi-R_curv.shape.gii
 smriprep/sub-01/anat/sub-01_hemi-R_inflated.surf.gii
 smriprep/sub-01/anat/sub-01_hemi-R_midthickness.surf.gii
 smriprep/sub-01/anat/sub-01_hemi-R_pial.surf.gii
+smriprep/sub-01/anat/sub-01_hemi-R_space-fsLR_den-32k_midthickness.surf.gii
+smriprep/sub-01/anat/sub-01_hemi-R_space-fsLR_den-32k_pial.surf.gii
+smriprep/sub-01/anat/sub-01_hemi-R_space-fsLR_den-32k_white.surf.gii
 smriprep/sub-01/anat/sub-01_hemi-R_space-fsLR_desc-msmsulc_sphere.surf.gii
 smriprep/sub-01/anat/sub-01_hemi-R_space-fsLR_desc-reg_sphere.surf.gii
 smriprep/sub-01/anat/sub-01_hemi-R_space-fsaverage_desc-reg_sphere.surf.gii

src/smriprep/workflows/anatomical.py

@@ -84,7 +84,7 @@
     init_morph_grayords_wf,
     init_msm_sulc_wf,
     init_refinement_wf,
-    init_resample_midthickness_wf,
+    init_resample_surfaces_wf,
     init_surface_derivatives_wf,
     init_surface_recon_wf,
 )
@@ -379,11 +379,23 @@ def init_anat_preproc_wf(
 
         if cifti_output:
             hcp_morphometrics_wf = init_hcp_morphometrics_wf(omp_nthreads=omp_nthreads)
-            resample_midthickness_wf = init_resample_midthickness_wf(grayord_density=cifti_output)
+            resample_surfaces_wf = init_resample_surfaces_wf(
+                surfaces=['white', 'pial', 'midthickness'],
+                grayord_density=cifti_output,
+            )
             morph_grayords_wf = init_morph_grayords_wf(
                 grayord_density=cifti_output, omp_nthreads=omp_nthreads
             )
 
+            ds_fsLR_surfaces_wf = init_ds_surfaces_wf(
+                output_dir=output_dir,
+                surfaces=['white', 'pial', 'midthickness'],
+                entities={
+                    'space': 'fsLR',
+                    'density': '32k' if cifti_output == '91k' else '59k',
+                },
+                name='ds_fsLR_surfaces_wf',
+            )
             ds_grayord_metrics_wf = init_ds_grayord_metrics_wf(
                 bids_root=bids_root,
                 output_dir=output_dir,
@@ -401,7 +413,9 @@ def init_anat_preproc_wf(
                 (surface_derivatives_wf, hcp_morphometrics_wf, [
                     ('outputnode.curv', 'inputnode.curv'),
                 ]),
-                (anat_fit_wf, resample_midthickness_wf, [
+                (anat_fit_wf, resample_surfaces_wf, [
+                    ('outputnode.white', 'inputnode.white'),
+                    ('outputnode.pial', 'inputnode.pial'),
                     ('outputnode.midthickness', 'inputnode.midthickness'),
                     (
                         f"outputnode.sphere_reg_{'msm' if msm_sulc else 'fsLR'}",
@@ -421,12 +435,20 @@ def init_anat_preproc_wf(
                     ('outputnode.thickness', 'inputnode.thickness'),
                     ('outputnode.roi', 'inputnode.roi'),
                 ]),
-                (resample_midthickness_wf, morph_grayords_wf, [
+                (resample_surfaces_wf, morph_grayords_wf, [
                     ('outputnode.midthickness_fsLR', 'inputnode.midthickness_fsLR'),
                 ]),
+                (anat_fit_wf, ds_fsLR_surfaces_wf, [
+                    ('outputnode.t1w_valid_list', 'inputnode.source_files'),
+                ]),
                 (anat_fit_wf, ds_grayord_metrics_wf, [
                     ('outputnode.t1w_valid_list', 'inputnode.source_files'),
                 ]),
+                (resample_surfaces_wf, ds_fsLR_surfaces_wf, [
+                    ('outputnode.white_fsLR', 'inputnode.white'),
+                    ('outputnode.pial_fsLR', 'inputnode.pial'),
+                    ('outputnode.midthickness_fsLR', 'inputnode.midthickness'),
+                ]),
                 (morph_grayords_wf, ds_grayord_metrics_wf, [
                     ('outputnode.curv_fsLR', 'inputnode.curv'),
                     ('outputnode.curv_metadata', 'inputnode.curv_metadata'),

src/smriprep/workflows/outputs.py

@@ -708,6 +708,7 @@ def init_ds_surfaces_wf(
     *,
     output_dir: str,
     surfaces: list[str],
+    entities: dict[str, str] | None = None,
     name='ds_surfaces_wf',
 ) -> Workflow:
     """
@@ -721,6 +722,8 @@ def init_ds_surfaces_wf(
         Directory in which to save derivatives
     surfaces : :class:`str`
         List of surfaces to generate DataSinks for
+    entities : :class:`dict` of :class:`str`
+        Entities to include in outputs
     name : :class:`str`
         Workflow name (default: ds_surfaces_wf)
 
@@ -739,6 +742,9 @@ def init_ds_surfaces_wf(
     """
     workflow = Workflow(name=name)
 
+    if entities is None:
+        entities = {}
+
     inputnode = pe.Node(
         niu.IdentityInterface(fields=['source_files'] + surfaces),
         name='inputnode',
@@ -766,6 +772,8 @@ def init_ds_surfaces_wf(
             elif surf == 'sphere_reg_msm':
                 ds_surf.inputs.space, ds_surf.inputs.desc = 'fsLR', 'msmsulc'
 
+        ds_surf.inputs.trait_set(**entities)
+
         # fmt:off
         workflow.connect([
             (inputnode, ds_surf, [(surf, 'in_file'), ('source_files', 'source_file')]),

src/smriprep/workflows/surfaces.py

@@ -1316,39 +1316,45 @@ def init_anat_ribbon_wf(name='anat_ribbon_wf'):
     return workflow
 
 
-def init_resample_midthickness_wf(
+def init_resample_surfaces_wf(
+    surfaces: list[str],
     grayord_density: ty.Literal['91k', '170k'],
-    name: str = 'resample_midthickness_wf',
+    name: str = 'resample_surfaces_wf',
 ):
     """
-    Resample subject midthickness surface to specified density.
+    Resample subject surfaces surface to specified density.
 
     Workflow Graph
         .. workflow::
             :graph2use: colored
             :simple_form: yes
 
-            from smriprep.workflows.surfaces import init_resample_midthickness_wf
-            wf = init_resample_midthickness_wf(grayord_density="91k")
+            from smriprep.workflows.surfaces import init_resample_surfaces_wf
+            wf = init_resample_surfaces_wf(
+                surfaces=['white', 'pial', 'midthickness'],
+                grayord_density='91k',
+            )
 
     Parameters
     ----------
-    grayord_density : :obj:`str`
+    surfaces : :class:`list` of :class:`str`
+        Names of surfaces (e.g., ``'white'``) to resample. Both hemispheres will be resampled.
+    grayord_density : :class:`str`
         Either `91k` or `170k`, representing the total of vertices or *grayordinates*.
-    name : :obj:`str`
-        Unique name for the subworkflow (default: ``"resample_midthickness_wf"``)
+    name : :class:`str`
+        Unique name for the subworkflow (default: ``"resample_surfaces_wf"``)
 
     Inputs
     ------
-    midthickness
-        GIFTI surface mesh corresponding to the midthickness surface
+    ``<surface>``
+        Left and right GIFTIs for each surface name passed to ``surfaces``
     sphere_reg_fsLR
         GIFTI surface mesh corresponding to the subject's fsLR registration sphere
 
     Outputs
     -------
-    midthickness
-        GIFTI surface mesh corresponding to the midthickness surface, resampled to fsLR
+    ``<surface>``
+        Left and right GIFTI surface mesh corresponding to the input surface, resampled to fsLR
     """
     import templateflow.api as tf
     from niworkflows.engine.workflows import LiterateWorkflow as Workflow
@@ -1358,11 +1364,19 @@ def init_resample_midthickness_wf(
     fslr_density = '32k' if grayord_density == '91k' else '59k'
 
     inputnode = pe.Node(
-        niu.IdentityInterface(fields=['midthickness', 'sphere_reg_fsLR']),
+        niu.IdentityInterface(fields=[*surfaces, 'sphere_reg_fsLR']),
         name='inputnode',
     )
 
-    outputnode = pe.Node(niu.IdentityInterface(fields=['midthickness_fsLR']), name='outputnode')
+    outputnode = pe.Node(
+        niu.IdentityInterface(fields=[f'{surf}_fsLR' for surf in surfaces]), name='outputnode'
+    )
+
+    surface_list = pe.Node(
+        niu.Merge(len(surfaces), ravel_inputs=True),
+        name='surface_list',
+        run_without_submitting=True,
+    )
 
     resampler = pe.MapNode(
         SurfaceResample(method='BARYCENTRIC'),
@@ -1380,15 +1394,30 @@ def init_resample_midthickness_wf(
                 extension='.surf.gii',
             )
         )
+        # Order matters. Iterate over surfaces, then hemis to get L R L R L R
+        for _surf in surfaces
         for hemi in ['L', 'R']
     ]
 
+    surface_groups = pe.Node(
+        niu.Split(splits=[2] * len(surfaces)),
+        name='surface_groups',
+        run_without_submitting=True,
+    )
+
     workflow.connect([
+        (inputnode, surface_list, [
+            ((surf, _sorted_by_basename), f'in{i}')
+            for i, surf in enumerate(surfaces, start=1)
+        ]),
         (inputnode, resampler, [
-            ('midthickness', 'surface_in'),
-            ('sphere_reg_fsLR', 'current_sphere'),
+            (('sphere_reg_fsLR', _repeat, len(surfaces)), 'current_sphere'),
+        ]),
+        (surface_list, resampler, [('out', 'surface_in')]),
+        (resampler, surface_groups, [('surface_out', 'inlist')]),
+        (surface_groups, outputnode, [
+            (f'out{i}', f'{surf}_fsLR') for i, surf in enumerate(surfaces, start=1)
         ]),
-        (resampler, outputnode, [('surface_out', 'midthickness_fsLR')]),
     ])  # fmt:skip
 
     return workflow
@@ -1678,3 +1707,7 @@ def _select_seg(in_files, segmentation):
             if segmentation in fl:
                 return fl
     raise FileNotFoundError(f'No segmentation containing "{segmentation}" was found.')
+
+
+def _repeat(seq: list, count: int) -> list:
+    return seq * count