Merge pull request #482 from tsalo/cortex-mask

ENH: Generate cortex mask in anat_fit_wf

Author: effigies

.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,32 @@ 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()
+        ds_cortex_masks_wf = init_ds_surface_masks_wf(
+            output_dir=output_dir,
+            mask_type='cortex',
+            name='ds_cortex_masks_wf',
+        )
+
+        workflow.connect([
+            (surfaces_buffer, cortex_masks_wf, [
+                ('midthickness', 'inputnode.midthickness'),
+                ('thickness', 'inputnode.thickness'),
+            ]),
+            (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,101 @@ 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.JoinNode(
+        niu.IdentityInterface(fields=['mask_files']), name='outputnode', joinsource='ds_itersource'
+    )
+
+    ds_itersource = pe.Node(
+        niu.IdentityInterface(fields=['hemi']),
+        name='ds_itersource',
+        iterables=[('hemi', ['L', 'R'])],
+    )
+
+    sources = pe.Node(niu.Function(function=_bids_relative), name='sources')
+    sources.inputs.bids_root = output_dir
+
+    select_files = pe.Node(
+        KeySelect(fields=['mask_file', 'sources'], keys=['L', 'R']),
+        name='select_files',
+        run_without_submitting=True,
+    )
+
+    ds_surf_mask = pe.Node(
+        DerivativesDataSink(
+            base_directory=output_dir,
+            suffix='mask',
+            desc=mask_type,
+            extension='.label.gii',
+            Type='Brain' if mask_type == 'brain' else 'ROI',
+            **entities,
+        ),
+        name='ds_surf_mask',
+        run_without_submitting=True,
+    )
+
+    workflow.connect([
+        (inputnode, select_files, [
+            ('mask_files', 'mask_file'),
+            ('source_files', 'sources'),
+        ]),
+        (select_files, sources, [('sources', 'in_files')]),
+        (ds_itersource, select_files, [('hemi', 'key')]),
+        (ds_itersource, ds_surf_mask, [('hemi', 'hemi')]),
+        (select_files, ds_surf_mask, [
+            ('mask_file', 'in_file'),
+            (('sources', _pop), 'source_file'),
+        ]),
+        (sources, ds_surf_mask, [('out', 'Sources')]),
+        (ds_surf_mask, outputnode, [('out_file', 'mask_files')]),
+    ])  # fmt: skip
+
+    return workflow
+
+
 def _bids_relative(in_files, bids_root):
     from pathlib import Path
 
@@ -1335,3 +1430,7 @@ def _read_json(in_file):
     from pathlib import Path
 
     return loads(Path(in_file).read_text())
+
+
+def _pop(in_list):
+    return in_list[0]

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,15 +1151,99 @@ 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')]),
+    ])  # 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',
+    )
+
+    itersource = pe.Node(
+        niu.IdentityInterface(fields=['hemi']),
+        name='itersource',
+        iterables=[('hemi', ['L', 'R'])],
+    )
+
+    outputnode = pe.JoinNode(
+        niu.IdentityInterface(fields=['cortex_masks', 'source_files']),
+        name='outputnode',
+        joinsource='itersource',
+    )
+
+    select_surfaces = pe.Node(
+        KeySelect(fields=['thickness', 'midthickness'], keys=['L', 'R']),
+        name='select_surfaces',
+        run_without_submitting=True,
+    )
+
+    combine_sources = pe.Node(niu.Merge(2), name='combine_sources', run_without_submitting=True)
+
+    abs_thickness = pe.Node(
+        MetricMath(metric='thickness', operation='abs'),
+        name='abs_thickness',
+        mem_gb=DEFAULT_MEMORY_MIN_GB,
+    )
+    initial_roi = pe.Node(
+        MetricMath(metric='roi', operation='bin'), name='initial_roi', mem_gb=DEFAULT_MEMORY_MIN_GB
+    )
+    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, select_surfaces, [
+            ('thickness', 'thickness'),
+            ('midthickness', 'midthickness'),
+        ]),
+        (itersource, select_surfaces, [('hemi', 'key')]),
+        (itersource, abs_thickness, [('hemi', 'hemisphere')]),
         (itersource, initial_roi, [('hemi', 'hemisphere')]),
-        (abs_thickness, initial_roi, [('metric_file', 'metric_file')]),
+        (select_surfaces, abs_thickness, [('thickness', 'metric_file')]),
         (select_surfaces, fill_holes, [('midthickness', 'surface_file')]),
         (select_surfaces, native_roi, [('midthickness', 'surface_file')]),
+        (abs_thickness, initial_roi, [('metric_file', 'metric_file')]),
         (initial_roi, fill_holes, [('metric_file', 'metric_file')]),
         (fill_holes, native_roi, [('out_file', 'metric_file')]),
-        (native_roi, outputnode, [('out_file', 'roi')]),
+        (native_roi, outputnode, [('out_file', 'cortex_masks')]),
+        (select_surfaces, combine_sources, [
+            ('midthickness', 'in1'),
+            ('thickness', 'in2'),
+        ]),
+        (combine_sources, outputnode, [('out', 'source_files')]),
     ])  # fmt:skip
 
     return workflow