diff --git a/.circleci/ds005_outputs.txt b/.circleci/ds005_outputs.txt index c96a92743c..2167c1bcc1 100644 --- a/.circleci/ds005_outputs.txt +++ b/.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 diff --git a/pyproject.toml b/pyproject.toml index cf8e0a5d6d..9f719c98f3 100644 --- a/pyproject.toml +++ b/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", diff --git a/src/smriprep/data/io_spec.json b/src/smriprep/data/io_spec.json index 54edaf3ac1..b3c97c588b 100644 --- a/src/smriprep/data/io_spec.json +++ b/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": { diff --git a/src/smriprep/workflows/anatomical.py b/src/smriprep/workflows/anatomical.py index 31107912df..2c1c3a9806 100644 --- a/src/smriprep/workflows/anatomical.py +++ b/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 diff --git a/src/smriprep/workflows/outputs.py b/src/smriprep/workflows/outputs.py index a2f946356e..9a750f145e 100644 --- a/src/smriprep/workflows/outputs.py +++ b/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] diff --git a/src/smriprep/workflows/surfaces.py b/src/smriprep/workflows/surfaces.py index 2840a166b6..3dce50fc81 100644 --- a/src/smriprep/workflows/surfaces.py +++ b/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