ENH: Generate cortex mask in anat_fit_wf

.circleci/ds005_outputs.txt

@@ -22,6 +22,8 @@ smriprep/sub-01/anat/sub-01_dseg.nii.gz
 smriprep/sub-01/anat/sub-01_from-fsnative_to-T1w_mode-image_xfm.txt
 smriprep/sub-01/anat/sub-01_from-T1w_to-fsnative_mode-image_xfm.txt
 smriprep/sub-01/anat/sub-01_hemi-L_curv.shape.gii
+smriprep/sub-01/anat/sub-01_hemi-L_desc-cortex_mask.json
+smriprep/sub-01/anat/sub-01_hemi-L_desc-cortex_mask.label.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
@@ -36,6 +38,8 @@ smriprep/sub-01/anat/sub-01_hemi-L_sulc.shape.gii
 smriprep/sub-01/anat/sub-01_hemi-L_thickness.shape.gii
 smriprep/sub-01/anat/sub-01_hemi-L_white.surf.gii
 smriprep/sub-01/anat/sub-01_hemi-R_curv.shape.gii
+smriprep/sub-01/anat/sub-01_hemi-R_desc-cortex_mask.json
+smriprep/sub-01/anat/sub-01_hemi-R_desc-cortex_mask.label.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

pyproject.toml

@@ -28,7 +28,7 @@ dependencies = [
     "nibabel >= 4.0.1",
     "nipype >= 1.8.5",
     "nireports >= 25.2.0",
-    "niworkflows >= 1.13.4",
+    "niworkflows @ git+https://github.com/nipreps/niworkflows.git@master",
     "numpy >= 1.24",
     "packaging >= 24",
     "pybids >= 0.16",

src/smriprep/data/io_spec.json

@@ -144,6 +144,14 @@
         "desc": "msmsulc",
         "suffix": "sphere",
         "extension": ".surf.gii"
+      },
+      "cortex_mask": {
+        "datatype": "anat",
+        "hemi": ["L", "R"],
+        "space": null,
+        "desc": "cortex",
+        "suffix": "mask",
+        "extension": ".label.gii"
       }
     },
     "masks": {

src/smriprep/workflows/anatomical.py

@@ -69,6 +69,7 @@
     init_ds_fs_segs_wf,
     init_ds_grayord_metrics_wf,
     init_ds_mask_wf,
+    init_ds_surface_masks_wf,
     init_ds_surface_metrics_wf,
     init_ds_surfaces_wf,
     init_ds_template_registration_wf,
@@ -78,6 +79,7 @@
 )
 from .surfaces import (
     init_anat_ribbon_wf,
+    init_cortex_masks_wf,
     init_fsLR_reg_wf,
     init_gifti_morphometrics_wf,
     init_gifti_surfaces_wf,
@@ -429,12 +431,12 @@ def init_anat_preproc_wf(
                         f"outputnode.sphere_reg_{'msm' if msm_sulc else 'fsLR'}",
                         'inputnode.sphere_reg_fsLR',
                     ),
+                    ('outputnode.cortex_mask', 'inputnode.roi'),
                 ]),
                 (hcp_morphometrics_wf, morph_grayords_wf, [
                     ('outputnode.curv', 'inputnode.curv'),
                     ('outputnode.sulc', 'inputnode.sulc'),
                     ('outputnode.thickness', 'inputnode.thickness'),
-                    ('outputnode.roi', 'inputnode.roi'),
                 ]),
                 (resample_surfaces_wf, morph_grayords_wf, [
                     ('outputnode.midthickness_fsLR', 'inputnode.midthickness_fsLR'),
@@ -668,6 +670,7 @@ def init_anat_fit_wf(
                 'sphere_reg',
                 'sphere_reg_fsLR',
                 'sphere_reg_msm',
+                'cortex_mask',
                 'anat_ribbon',
                 # Reverse transform; not computable from forward transform
                 'std2anat_xfm',
@@ -1344,6 +1347,35 @@ def init_anat_fit_wf(
     else:
         LOGGER.info('ANAT Stage 10: MSM-Sulc disabled')
 
+    # Stage 11: Cortical surface mask
+    if len(precomputed.get('cortex_mask', [])) < 2:
+        LOGGER.info('ANAT Stage 11: Creating cortical surface mask')
+
+        cortex_masks_wf = init_cortex_masks_wf()
+        workflow.connect([
+            (surfaces_buffer, cortex_masks_wf, [
+                ('midthickness', 'inputnode.midthickness'),
+                ('thickness', 'inputnode.thickness'),
+            ]),
+        ])  # fmt:skip
+
+        ds_cortex_masks_wf = init_ds_surface_masks_wf(
+            output_dir=output_dir,
+            mask_type='cortex',
+            name='ds_cortex_masks_wf',
+            entities={'extension': '.label.gii'},
+        )
+        workflow.connect([
+            (cortex_masks_wf, ds_cortex_masks_wf, [
+                ('outputnode.cortex_masks', 'inputnode.mask_files'),
+                ('outputnode.source_files', 'inputnode.source_files'),
+            ]),
+            (ds_cortex_masks_wf, outputnode, [('outputnode.mask_files', 'cortex_mask')]),
+        ])  # fmt:skip
+    else:
+        LOGGER.info('ANAT Stage 11: Found pre-computed cortical surface mask')
+        outputnode.inputs.cortex_mask = sorted(precomputed['cortex_mask'])
+
     return workflow
 
 

src/smriprep/workflows/outputs.py

@@ -1231,6 +1231,105 @@ def init_template_iterator_wf(
     return workflow
 
 
+def init_ds_surface_masks_wf(
+    *,
+    output_dir: str,
+    mask_type: ty.Literal['cortex', 'roi', 'ribbon', 'brain'],
+    entities: dict[str, str] | None = None,
+    name='ds_surface_masks_wf',
+) -> Workflow:
+    """Save GIFTI surface masks.
+
+    Parameters
+    ----------
+    output_dir : :class:`str`
+        Directory in which to save derivatives
+    mask_type : :class:`str`
+        Type of mask to save
+    entities : :class:`dict` of :class:`str`
+        Entities to include in outputs
+    name : :class:`str`
+        Workflow name (default: ds_surface_masks_wf)
+
+    Inputs
+    ------
+    source_files : list of lists of str
+        List of lists of source files.
+        Left hemisphere sources first, then right hemisphere sources.
+    mask_files : list of str
+        List of input mask files.
+        Left hemisphere mask first, then right hemisphere mask.
+
+    Outputs
+    -------
+    mask_files : list of str
+        List of output mask files.
+        Left hemisphere mask first, then right hemisphere mask.
+    """
+    workflow = Workflow(name=name)
+
+    if entities is None:
+        entities = {}
+
+    inputnode = pe.Node(
+        niu.IdentityInterface(fields=['mask_files', 'source_files']),
+        name='inputnode',
+    )
+    outputnode = pe.Node(niu.IdentityInterface(fields=['mask_files']), name='outputnode')
+
+    combine_masks = pe.Node(
+        niu.Merge(2),
+        name='combine_masks',
+    )
+    workflow.connect([(combine_masks, outputnode, [('out', 'mask_files')])])
+
+    for i_hemi, hemi in enumerate(['L', 'R']):
+        select_mask = pe.Node(
+            niu.Select(index=i_hemi),
+            name=f'select_mask_{hemi}',
+            run_without_submitting=True,
+        )
+        workflow.connect([(inputnode, select_mask, [('mask_files', 'inlist')])])
+
+        select_source = pe.Node(
+            niu.Select(index=i_hemi),
+            name=f'select_source_{hemi}',
+            run_without_submitting=True,
+        )
+        workflow.connect([(inputnode, select_source, [('source_files', 'inlist')])])
+
+        sources = pe.Node(
+            niu.Function(function=_bids_relative),
+            name=f'sources_{hemi}',
+        )
+        sources.inputs.bids_root = output_dir
+
+        ds_mask = pe.Node(
+            DerivativesDataSink(
+                base_directory=output_dir,
+                hemi=hemi,
+                desc=mask_type,
+                **entities,
+            ),
+            name=f'ds_mask_{hemi}',
+            run_without_submitting=True,
+        )
+        if mask_type == 'brain':
+            ds_mask.inputs.Type = 'Brain'
+        else:
+            ds_mask.inputs.Type = 'ROI'
+
+        workflow.connect([
+            (select_mask, ds_mask, [('out', 'in_file')]),
+            (select_source, sources, [('out', 'in_files')]),
+            (select_source, ds_mask, [('out', 'source_file')]),
+            (sources, ds_mask, [('out', 'Sources')]),
+            (ds_mask, combine_masks, [('out_file', f'in{i_hemi + 1}')]),
+        ])  # fmt:skip
+
+    return workflow
+
+
 def _bids_relative(in_files, bids_root):
     from pathlib import Path
 

src/smriprep/workflows/surfaces.py

@@ -1071,8 +1071,6 @@ def init_hcp_morphometrics_wf(
         HCP-style curvature file in GIFTI format
     sulc
         HCP-style sulcal depth file in GIFTI format
-    roi
-        HCP-style cortical ROI file in GIFTI format
     """
     DEFAULT_MEMORY_MIN_GB = 0.01
 
@@ -1090,7 +1088,7 @@ def init_hcp_morphometrics_wf(
     )
 
     outputnode = pe.JoinNode(
-        niu.IdentityInterface(fields=['thickness', 'curv', 'sulc', 'roi']),
+        niu.IdentityInterface(fields=['thickness', 'curv', 'sulc']),
         name='outputnode',
         joinsource='itersource',
     )
@@ -1115,11 +1113,6 @@ def init_hcp_morphometrics_wf(
     # Thickness is presumably already positive, but HCP uses abs(-thickness)
     abs_thickness = pe.Node(MetricMath(metric='thickness', operation='abs'), name='abs_thickness')
 
-    # Native ROI is thickness > 0, with holes and islands filled
-    initial_roi = pe.Node(MetricMath(metric='roi', operation='bin'), name='initial_roi')
-    fill_holes = pe.Node(MetricFillHoles(), name='fill_holes', mem_gb=DEFAULT_MEMORY_MIN_GB)
-    native_roi = pe.Node(MetricRemoveIslands(), name='native_roi', mem_gb=DEFAULT_MEMORY_MIN_GB)
-
     # Dilation happens separately from ROI creation
     dilate_curv = pe.Node(
         MetricDilate(distance=10, nearest=True),
@@ -1158,20 +1151,120 @@ def init_hcp_morphometrics_wf(
         (dilate_curv, outputnode, [('out_file', 'curv')]),
         (dilate_thickness, outputnode, [('out_file', 'thickness')]),
         (invert_sulc, outputnode, [('metric_file', 'sulc')]),
-        # Native ROI file from thickness
-        (inputnode, initial_roi, [('subject_id', 'subject_id')]),
-        (itersource, initial_roi, [('hemi', 'hemisphere')]),
-        (abs_thickness, initial_roi, [('metric_file', 'metric_file')]),
-        (select_surfaces, fill_holes, [('midthickness', 'surface_file')]),
-        (select_surfaces, native_roi, [('midthickness', 'surface_file')]),
-        (initial_roi, fill_holes, [('metric_file', 'metric_file')]),
-        (fill_holes, native_roi, [('out_file', 'metric_file')]),
-        (native_roi, outputnode, [('out_file', 'roi')]),
     ])  # fmt:skip
 
     return workflow
 
 
+def init_cortex_masks_wf(
+    *,
+    name: str = 'cortex_masks_wf',
+):
+    """Create cortical surface masks from surface files.
+
+    Workflow Graph
+        .. workflow::
+            :graph2use: orig
+            :simple_form: yes
+
+            from smriprep.workflows.surfaces import init_cortex_masks_wf
+            wf = init_cortex_masks_wf()
+
+    Inputs
+    ------
+    midthickness : len-2 list of str
+        Each hemisphere's FreeSurfer midthickness surface file in GIFTI format
+    thickness : len-2 list of str
+        Each hemisphere's FreeSurfer thickness file in GIFTI format
+
+    Outputs
+    -------
+    cortex_masks : len-2 list of str
+        Cortical surface mask in GIFTI format for each hemisphere
+    source_files : len-2 list of lists of str
+        Each hemisphere's source files, which are used to create the mask
+    """
+    DEFAULT_MEMORY_MIN_GB = 0.01
+
+    workflow = Workflow(name=name)
+
+    inputnode = pe.Node(
+        niu.IdentityInterface(fields=['midthickness', 'thickness']),
+        name='inputnode',
+    )
+    outputnode = pe.Node(
+        niu.IdentityInterface(fields=['cortex_masks', 'source_files']),
+        name='outputnode',
+    )
+
+    # Combine the inputs into a list
+    combine_sources = pe.Node(
+        niu.Merge(2, no_flatten=True),
+        name='combine_sources',
+    )
+    workflow.connect([
+        (inputnode, combine_sources, [
+            ('midthickness', 'in1'),
+            ('thickness', 'in2'),
+        ]),
+        (combine_sources, outputnode, [(('out', _transpose_lol), 'source_files')]),
+    ])  # fmt:skip
+
+    combine_masks = pe.Node(
+        niu.Merge(2),
+        name='combine_masks',
+    )
+    workflow.connect([(combine_masks, outputnode, [('out', 'cortex_masks')])])
+
+    for i_hemi, hemi in enumerate(['L', 'R']):
+        select_midthickness = pe.Node(
+            niu.Select(index=i_hemi),
+            name=f'select_midthickness_{hemi}',
+        )
+        select_thickness = pe.Node(
+            niu.Select(index=i_hemi),
+            name=f'select_thickness_{hemi}',
+        )
+        workflow.connect([
+            (inputnode, select_midthickness, [('midthickness', 'inlist')]),
+            (inputnode, select_thickness, [('thickness', 'inlist')]),
+        ])  # fmt:skip
+
+        # Thickness is presumably already positive, but HCP uses abs(-thickness)
+        abs_thickness = pe.Node(
+            MetricMath(metric='thickness', operation='abs', hemisphere=hemi),
+            name=f'abs_thickness_{hemi}',
+        )
+
+        # Native ROI is thickness > 0, with holes and islands filled
+        initial_roi = pe.Node(
+            MetricMath(metric='roi', operation='bin', hemisphere=hemi),
+            name=f'initial_roi_{hemi}',
+        )
+        fill_holes = pe.Node(
+            MetricFillHoles(),
+            name=f'fill_holes_{hemi}',
+            mem_gb=DEFAULT_MEMORY_MIN_GB,
+        )
+        native_roi = pe.Node(
+            MetricRemoveIslands(),
+            name=f'native_roi_{hemi}',
+            mem_gb=DEFAULT_MEMORY_MIN_GB,
+        )
+
+        workflow.connect([
+            (select_thickness, abs_thickness, [('out', 'metric_file')]),
+            (abs_thickness, initial_roi, [('metric_file', 'metric_file')]),
+            (select_midthickness, fill_holes, [('out', 'surface_file')]),
+            (select_midthickness, native_roi, [('out', 'surface_file')]),
+            (initial_roi, fill_holes, [('metric_file', 'metric_file')]),
+            (fill_holes, native_roi, [('out_file', 'metric_file')]),
+            (native_roi, combine_masks, [('out_file', f'in{i_hemi + 1}')]),
+        ])  # fmt:skip
+
+    return workflow
+
+
 @tag('anat.segs-to-anat')
 def init_segs_to_native_wf(
     *,
@@ -1711,3 +1804,8 @@ def _select_seg(in_files, segmentation):
 
 def _repeat(seq: list, count: int) -> list:
     return seq * count
+
+
+def _transpose_lol(inlist):
+    """Transpose a list of lists."""
+    return list(map(list, zip(*inlist, strict=False)))