ENH: Generate cortex mask in anat_fit_wf

.circleci/ds005_outputs.txt

@@ -22,6 +22,7 @@ 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.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 +37,7 @@ 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.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

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

@@ -78,6 +78,7 @@
 )
 from .surfaces import (
     init_anat_ribbon_wf,
+    init_cortex_mask_wf,
     init_fsLR_reg_wf,
     init_gifti_morphometrics_wf,
     init_gifti_surfaces_wf,
@@ -429,12 +430,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 +669,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 +1346,65 @@ 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')
+
+        # Merge outputs into a single list
+        merge_cortex_masks = pe.Node(
+            niu.Merge(2),
+            name='merge_cortex_masks',
+        )
+        workflow.connect([(merge_cortex_masks, outputnode, [('out', 'cortex_mask')])])
+
+        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}',
+            )
+            cortex_mask_wf = init_cortex_mask_wf(name=f'cortex_mask_wf_{hemi}')
+            cortex_mask_wf.inputs.inputnode.hemi = hemi
+
+            workflow.connect([
+                (surfaces_buffer, select_midthickness, [('midthickness', 'inlist')]),
+                (surfaces_buffer, select_thickness, [('thickness', 'inlist')]),
+                (select_midthickness, cortex_mask_wf, [('out', 'inputnode.midthickness')]),
+                (select_thickness, cortex_mask_wf, [('out', 'inputnode.thickness')]),
+            ])  # fmt:skip
+
+            # Combine the inputs into a list
+            combine_inputs = pe.Node(
+                niu.Merge(2),
+                name=f'combine_inputs_{hemi}',
+            )
+            workflow.connect([
+                (select_midthickness, combine_inputs, [('out', 'in1')]),
+                (select_thickness, combine_inputs, [('out', 'in2')]),
+            ])  # fmt:skip
+
+            ds_cortex_mask_wf = init_ds_mask_wf(
+                bids_root=bids_root,
+                output_dir=output_dir,
+                mask_type='roi',
+                name=f'ds_cortex_mask_wf_{hemi}',
+                extra_entities={'extension': '.label.gii', 'hemi': hemi},
+            )
+
+            workflow.connect([
+                (cortex_mask_wf, ds_cortex_mask_wf, [('outputnode.roi', 'inputnode.mask_file')]),
+                (combine_inputs, ds_cortex_mask_wf, [('out', 'inputnode.source_files')]),
+                (ds_cortex_mask_wf, merge_cortex_masks, [
+                    ('outputnode.mask_file', f'in{i_hemi + 1}'),
+                ]),
+            ])  # 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/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,12 +1151,66 @@ 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')]),
+    ])  # fmt:skip
+
+    return workflow
+
+
+def init_cortex_mask_wf(
+    *,
+    name: str = 'cortex_mask_wf',
+):
+    """Create a cortical surface mask from a surface file.
+
+    Workflow Graph
+        .. workflow::
+            :graph2use: orig
+            :simple_form: yes
+
+            from smriprep.workflows.surfaces import init_cortex_mask_wf
+            wf = init_cortex_mask_wf()
+
+    Inputs
+    ------
+    midthickness : str
+        One hemisphere's FreeSurfer midthickness surface file in GIFTI format
+    thickness : str
+        One hemisphere's FreeSurfer thickness file in GIFTI format
+    hemi : {'L', 'R'}
+        Hemisphere indicator
+
+    Outputs
+    -------
+    roi : str
+        Cortical surface mask in GIFTI format
+    """
+    DEFAULT_MEMORY_MIN_GB = 0.01
+
+    workflow = Workflow(name=name)
+
+    inputnode = pe.Node(
+        niu.IdentityInterface(fields=['midthickness', 'thickness', 'hemi']),
+        name='inputnode',
+    )
+    outputnode = pe.Node(niu.IdentityInterface(fields=['roi']), name='outputnode')
+
+    # 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)
+
+    workflow.connect([
+        (inputnode, abs_thickness, [
+            ('hemi', 'hemisphere'),
+            ('thickness', 'metric_file'),
+        ]),
+        (inputnode, 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')]),
+        (inputnode, fill_holes, [('midthickness', 'surface_file')]),
+        (inputnode, 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')]),