diff --git a/.circleci/config.yml b/.circleci/config.yml index c45ce9b242..002a6a2af8 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -12,6 +12,9 @@ _machine_defaults: &machine_defaults _python_defaults: &python_defaults docker: - image: cimg/python:3.10.9 + auth: + username: $DOCKER_USER + password: $DOCKER_PAT working_directory: /tmp/src/smriprep _docker_auth: &docker_auth @@ -43,6 +46,9 @@ _pull_from_registry: &pull_from_registry docs_deploy: &docs docker: - image: node:8.10.0 + auth: + username: $DOCKER_USER + password: $DOCKER_PAT working_directory: /tmp/gh-pages steps: - run: diff --git a/.circleci/ds005_outputs.txt b/.circleci/ds005_outputs.txt index 3e26840d02..779c44abc6 100644 --- a/.circleci/ds005_outputs.txt +++ b/.circleci/ds005_outputs.txt @@ -16,6 +16,7 @@ smriprep/sub-01/anat/sub-01_desc-brain_mask.json smriprep/sub-01/anat/sub-01_desc-brain_mask.nii.gz smriprep/sub-01/anat/sub-01_desc-preproc_T1w.json smriprep/sub-01/anat/sub-01_desc-preproc_T1w.nii.gz +smriprep/sub-01/anat/sub-01_desc-ribbon_mask.json smriprep/sub-01/anat/sub-01_desc-ribbon_mask.nii.gz smriprep/sub-01/anat/sub-01_dseg.nii.gz smriprep/sub-01/anat/sub-01_from-fsnative_to-T1w_mode-image_xfm.txt @@ -27,6 +28,7 @@ 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_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_sphere.surf.gii 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 @@ -37,6 +39,7 @@ 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_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_sphere.surf.gii smriprep/sub-01/anat/sub-01_hemi-R_sulc.shape.gii smriprep/sub-01/anat/sub-01_hemi-R_thickness.shape.gii smriprep/sub-01/anat/sub-01_hemi-R_white.surf.gii diff --git a/smriprep/data/io_spec.json b/smriprep/data/io_spec.json index b77d61276d..7a6382b882 100644 --- a/smriprep/data/io_spec.json +++ b/smriprep/data/io_spec.json @@ -70,68 +70,91 @@ "mode": "image" } }, - "std_xfms": { - "anat2std_xfm": { + "surfaces": { + "white": { "datatype": "anat", - "extension": "h5", - "from": "T1w", - "to": [], - "suffix": "xfm", - "mode": "image" + "hemi": ["L", "R"], + "space": null, + "suffix": "white", + "extension": ".surf.gii" }, - "std2anat_xfm": { + "pial": { "datatype": "anat", - "extension": "h5", - "from": [], - "to": "T1w", - "suffix": "xfm", - "mode": "image" - } - }, - "surfaces": { - "t1w_aseg": { + "hemi": ["L", "R"], + "space": null, + "suffix": "pial", + "extension": ".surf.gii" + }, + "midthickness": { "datatype": "anat", - "desc": "aseg", - "suffix": "dseg" + "hemi": ["L", "R"], + "space": null, + "suffix": "midthickness", + "extension": ".surf.gii" }, - "t1w_aparc": { + "sphere": { "datatype": "anat", - "desc": "aparcaseg", - "suffix": "dseg" + "hemi": ["L", "R"], + "space": null, + "desc": null, + "suffix": "sphere", + "extension": ".surf.gii" }, - "t1w2fsnative_xfm": { + "thickness": { "datatype": "anat", - "from": "T1w", - "to": "fsnative", - "extension": "txt", - "suffix": "xfm", - "mode": "image" + "hemi": ["L", "R"], + "space": null, + "suffix": "thickness", + "extension": ".shape.gii" }, - "fsnative2t1w_xfm": { + "sulc": { "datatype": "anat", - "from": "fsnative", - "to": "T1w", - "extension": "txt", - "suffix": "xfm", - "mode": "image" + "hemi": ["L", "R"], + "space": null, + "suffix": "sulc", + "extension": ".shape.gii" }, - "surfaces": { + "curv": { "datatype": "anat", "hemi": ["L", "R"], - "extension": "surf.gii", - "suffix": [ "inflated", "midthickness", "pial", "white"] + "space": null, + "suffix": "curv", + "extension": ".shape.gii" }, - "morphometrics": { + "sphere_reg": { "datatype": "anat", "hemi": ["L", "R"], - "extension": "shape.gii", - "suffix": ["thickness", "sulc", "curv"] + "space": null, + "desc": "reg", + "suffix": "sphere", + "extension": ".surf.gii" }, + "sphere_reg_fsLR": { + "datatype": "anat", + "hemi": ["L", "R"], + "space": "fsLR", + "desc": "reg", + "suffix": "sphere", + "extension": ".surf.gii" + }, + "sphere_reg_msm": { + "datatype": "anat", + "hemi": ["L", "R"], + "space": "fsLR", + "desc": "msmsulc", + "suffix": "sphere", + "extension": ".surf.gii" + } + }, + "masks": { "anat_ribbon": { "datatype": "anat", "desc": "ribbon", "suffix": "mask", - "extension": "nii.gz" + "extension": [ + ".nii.gz", + ".nii" + ] } } }, diff --git a/smriprep/utils/bids.py b/smriprep/utils/bids.py index 125ae478a5..5b9174d6dd 100644 --- a/smriprep/utils/bids.py +++ b/smriprep/utils/bids.py @@ -21,39 +21,13 @@ # https://www.nipreps.org/community/licensing/ # """Utilities to handle BIDS inputs.""" -from collections import defaultdict from pathlib import Path from json import loads from pkg_resources import resource_filename as pkgrf from bids.layout import BIDSLayout -from bids.layout.writing import build_path -def get_outputnode_spec(): - """ - Generate outputnode's fields from I/O spec file. - - Examples - -------- - >>> get_outputnode_spec() # doctest: +NORMALIZE_WHITESPACE - ['t1w_preproc', 't1w_mask', 't1w_dseg', 't1w_tpms', - 'std_preproc', 'std_mask', 'std_dseg', 'std_tpms', - 'anat2std_xfm', 'std2anat_xfm', - 't1w_aseg', 't1w_aparc', - 't1w2fsnative_xfm', 'fsnative2t1w_xfm', - 'surfaces', 'morphometrics', 'anat_ribbon'] - - """ - spec = loads(Path(pkgrf("smriprep", "data/io_spec.json")).read_text())["queries"] - fields = ["_".join((m, s)) for m in ("t1w", "std") for s in spec["baseline"].keys()] - fields += [s for s in spec["std_xfms"].keys()] - fields += [s for s in spec["surfaces"].keys()] - return fields - - -def collect_derivatives( - derivatives_dir, subject_id, std_spaces, freesurfer, spec=None, patterns=None -): +def collect_derivatives(derivatives_dir, subject_id, std_spaces, spec=None, patterns=None): """Gather existing derivatives and compose a cache.""" if spec is None or patterns is None: _spec, _patterns = tuple( @@ -65,27 +39,9 @@ def collect_derivatives( if patterns is None: patterns = _patterns - derivs_cache = defaultdict(list, {}) layout = BIDSLayout(derivatives_dir, config=["bids", "derivatives"], validate=False) - derivatives_dir = Path(derivatives_dir) - - def _check_item(item): - if not item: - return None - - if isinstance(item, str): - item = [item] - - result = [] - for i in item: - if not (derivatives_dir / i).exists(): - i = i.rstrip(".gz") - if not (derivatives_dir / i).exists(): - return None - result.append(str(derivatives_dir / i)) - - return result + derivs_cache = {} for k, q in spec["baseline"].items(): q["subject"] = subject_id item = layout.get(return_type='filename', **q) @@ -94,18 +50,6 @@ def _check_item(item): derivs_cache["t1w_%s" % k] = item[0] if len(item) == 1 else item - for space in std_spaces: - for k, q in spec["std_xfms"].items(): - q["subject"] = subject_id - q["from"] = q["from"] or space - q["to"] = q["to"] or space - item = layout.get(return_type='filename', **q) - if not item: - continue - derivs_cache[k] += item - - derivs_cache = dict(derivs_cache) # Back to a standard dictionary - transforms = derivs_cache.setdefault('transforms', {}) for space in std_spaces: for k, q in spec["transforms"].items(): @@ -118,18 +62,14 @@ def _check_item(item): continue transforms.setdefault(space, {})[k] = item[0] if len(item) == 1 else item - if freesurfer: - for k, q in spec["surfaces"].items(): - q["subject"] = subject_id - item = _check_item(build_path(q, patterns)) - if not item: - continue + for k, q in spec["surfaces"].items(): + q["subject"] = subject_id + item = layout.get(return_type='filename', **q) + if not item or len(item) != 2: + continue - if len(item) == 1: - item = item[0] - derivs_cache[k] = item + derivs_cache[k] = sorted(item) - derivs_cache["template"] = std_spaces return derivs_cache diff --git a/smriprep/workflows/anatomical.py b/smriprep/workflows/anatomical.py index fba9c9b76a..23139ade94 100644 --- a/smriprep/workflows/anatomical.py +++ b/smriprep/workflows/anatomical.py @@ -55,6 +55,7 @@ from .fit.registration import init_register_template_wf from .outputs import ( init_anat_reports_wf, + init_ds_surface_metrics_wf, init_ds_surfaces_wf, init_ds_template_wf, init_ds_mask_wf, @@ -66,6 +67,10 @@ ) from .surfaces import ( init_anat_ribbon_wf, + init_fsLR_reg_wf, + init_gifti_morphometrics_wf, + init_gifti_surfaces_wf, + init_msm_sulc_wf, init_surface_derivatives_wf, init_surface_recon_wf, init_refinement_wf, @@ -244,6 +249,7 @@ def init_anat_preproc_wf( freesurfer=freesurfer, hires=hires, longitudinal=longitudinal, + msm_sulc=msm_sulc, skull_strip_mode=skull_strip_mode, skull_strip_template=skull_strip_template, spaces=spaces, @@ -283,6 +289,9 @@ def init_anat_preproc_wf( ("outputnode.anat2std_xfm", "anat2std_xfm"), ("outputnode.std2anat_xfm", "std2anat_xfm"), ("outputnode.fsnative2t1w_xfm", "fsnative2t1w_xfm"), + ("outputnode.sphere_reg", "sphere_reg"), + (f"outputnode.sphere_reg_{'msm' if msm_sulc else 'fsLR'}", "sphere_reg_fsLR"), + ("outputnode.anat_ribbon", "anat_ribbon"), ]), (anat_fit_wf, anat_second_derivatives_wf, [ ('outputnode.template', 'inputnode.template'), @@ -296,18 +305,15 @@ def init_anat_preproc_wf( ]) # fmt:on if freesurfer: - surfaces = ["white", "pial", "midthickness", "inflated", "sphere_reg", "sphere_reg_fsLR"] - if msm_sulc: - surfaces.append("sphere_reg_msm") surface_derivatives_wf = init_surface_derivatives_wf( - msm_sulc=msm_sulc, cifti_output=cifti_output, - sloppy=sloppy, ) - anat_ribbon_wf = init_anat_ribbon_wf() ds_surfaces_wf = init_ds_surfaces_wf( - bids_root=bids_root, output_dir=output_dir, surfaces=surfaces + bids_root=bids_root, output_dir=output_dir, surfaces=["inflated"] + ) + ds_curv_wf = init_ds_surface_metrics_wf( + bids_root=bids_root, output_dir=output_dir, metrics=["curv"], name="ds_curv_wf" ) # fmt:off @@ -318,22 +324,19 @@ def init_anat_preproc_wf( ('outputnode.subject_id', 'inputnode.subject_id'), ('outputnode.fsnative2t1w_xfm', 'inputnode.fsnative2t1w_xfm'), ]), - (anat_fit_wf, anat_ribbon_wf, [ - ('outputnode.t1w_mask', 'inputnode.ref_file'), - ]), - (surface_derivatives_wf, anat_ribbon_wf, [ - ('outputnode.white', 'inputnode.white'), - ('outputnode.pial', 'inputnode.pial'), - ]), (anat_fit_wf, ds_surfaces_wf, [ ('outputnode.t1w_valid_list', 'inputnode.source_files'), ]), (surface_derivatives_wf, ds_surfaces_wf, [ - (f'outputnode.{surface}', f'inputnode.{surface}') - for surface in surfaces + ('outputnode.inflated', 'inputnode.inflated'), + ]), + (anat_fit_wf, ds_curv_wf, [ + ('outputnode.t1w_valid_list', 'inputnode.source_files'), + ]), + (surface_derivatives_wf, ds_curv_wf, [ + ('outputnode.curv', 'inputnode.curv'), ]), (surface_derivatives_wf, anat_second_derivatives_wf, [ - ('outputnode.morphometrics', 'inputnode.morphometrics'), ('outputnode.out_aseg', 'inputnode.t1w_fs_aseg'), ('outputnode.out_aparc', 'inputnode.t1w_fs_aparc'), ('outputnode.cifti_morph', 'inputnode.cifti_morph'), @@ -342,14 +345,6 @@ def init_anat_preproc_wf( (surface_derivatives_wf, outputnode, [ ('outputnode.out_aseg', 't1w_aseg'), ('outputnode.out_aparc', 't1w_aparc'), - ('outputnode.sphere_reg', 'sphere_reg'), - ('outputnode.sphere_reg_fsLR', 'sphere_reg_fsLR'), - ]), - (anat_ribbon_wf, outputnode, [ - ("outputnode.anat_ribbon", "anat_ribbon"), - ]), - (anat_ribbon_wf, anat_second_derivatives_wf, [ - ("outputnode.anat_ribbon", "inputnode.anat_ribbon"), ]), ]) # fmt:on @@ -364,6 +359,7 @@ def init_anat_fit_wf( freesurfer: bool, hires: bool, longitudinal: bool, + msm_sulc: bool, t1w: list, t2w: list, skull_strip_mode: str, @@ -402,6 +398,7 @@ def init_anat_fit_wf( freesurfer=True, hires=True, longitudinal=False, + msm_sulc=True, t1w=['t1w.nii.gz'], t2w=['t2w.nii.gz'], skull_strip_mode='force', @@ -546,6 +543,17 @@ def init_anat_fit_wf( "t1w_tpms", "anat2std_xfm", "fsnative2t1w_xfm", + # Surface and metric derivatives for fsLR resampling + "white", + "pial", + "midthickness", + "sphere", + "thickness", + "sulc", + "sphere_reg", + "sphere_reg_fsLR", + "sphere_reg_msm", + "anat_ribbon", # Reverse transform; not computable from forward transform "std2anat_xfm", # Metadata @@ -571,11 +579,6 @@ def init_anat_fit_wf( niu.IdentityInterface(fields=["t1w_preproc", "t1w_mask", "t1w_brain", "ants_seg"]), name="t1w_buffer", ) - # Refined stage 2 results; may be direct copy if no refinement - refined_buffer = pe.Node( - niu.IdentityInterface(fields=["t1w_mask", "t1w_brain"]), - name="refined_buffer", - ) # Stage 3 results seg_buffer = pe.Node( niu.IdentityInterface(fields=["t1w_dseg", "t1w_tpms"]), @@ -586,6 +589,24 @@ def init_anat_fit_wf( anat2std_buffer = pe.Node(niu.Merge(2), name="anat2std_buffer") std2anat_buffer = pe.Node(niu.Merge(2), name="std2anat_buffer") + # Stage 6 results: Refined stage 2 results; may be direct copy if no refinement + refined_buffer = pe.Node( + niu.IdentityInterface(fields=["t1w_mask", "t1w_brain"]), + name="refined_buffer", + ) + + # Stage 8 results: GIFTI surfaces + surfaces_buffer = pe.Node( + niu.IdentityInterface( + fields=["white", "pial", "midthickness", "sphere", "sphere_reg", "thickness", "sulc"] + ), + name="surfaces_buffer", + ) + + # Stage 9 and 10 results: fsLR sphere registration + fsLR_buffer = pe.Node(niu.IdentityInterface(fields=["sphere_reg_fsLR"]), name="fsLR_buffer") + msm_buffer = pe.Node(niu.IdentityInterface(fields=["sphere_reg_msm"]), name="msm_buffer") + # fmt:off workflow.connect([ (seg_buffer, outputnode, [ @@ -596,6 +617,17 @@ def init_anat_fit_wf( (std2anat_buffer, outputnode, [("out", "std2anat_xfm")]), (template_buffer, outputnode, [("out", "template")]), (sourcefile_buffer, outputnode, [("source_files", "t1w_valid_list")]), + (surfaces_buffer, outputnode, [ + ("white", "white"), + ("pial", "pial"), + ("midthickness", "midthickness"), + ("sphere", "sphere"), + ("sphere_reg", "sphere_reg"), + ("thickness", "thickness"), + ("sulc", "sulc"), + ]), + (fsLR_buffer, outputnode, [("sphere_reg_fsLR", "sphere_reg_fsLR")]), + (msm_buffer, outputnode, [("sphere_reg_msm", "sphere_reg_msm")]), ]) # fmt:on @@ -624,7 +656,7 @@ def init_anat_fit_wf( # If desc-preproc_T1w.nii.gz is provided, just validate it anat_validate = pe.Node(ValidateImage(), name="anat_validate", run_without_submitting=True) if not have_t1w: - LOGGER.info("Stage 1: Adding template workflow") + LOGGER.info("ANAT Stage 1: Adding template workflow") ants_ver = ANTsInfo.version() or "(version unknown)" desc += f"""\ {"Each" if num_t1w > 1 else "The"} T1w image was corrected for intensity @@ -660,7 +692,7 @@ def init_anat_fit_wf( ]) # fmt:on else: - LOGGER.info("Found preprocessed T1w - skipping Stage 1") + LOGGER.info("ANAT Found preprocessed T1w - skipping Stage 1") desc += """ A preprocessed T1w image was provided as a precomputed input and used as T1w-reference throughout the workflow. """ @@ -679,7 +711,7 @@ def init_anat_fit_wf( # We always need to generate t1w_brain; how to do that depends on whether we have # a pre-corrected T1w or precomputed mask, or are given an already masked image if not have_mask: - LOGGER.info("Stage 2: Preparing brain extraction workflow") + LOGGER.info("ANAT Stage 2: Preparing brain extraction workflow") if skull_strip_mode == "auto": run_skull_strip = all(_is_skull_stripped(img) for img in t1w) else: @@ -717,7 +749,7 @@ def init_anat_fit_wf( # fmt:on # Determine mask from T1w and uniformize elif not have_t1w: - LOGGER.info("Stage 2: Skipping skull-strip, INU-correction only") + LOGGER.info("ANAT Stage 2: Skipping skull-strip, INU-correction only") desc += """\ The provided T1w image was previously skull-stripped; a brain mask was derived from the input image. @@ -739,7 +771,7 @@ def init_anat_fit_wf( # fmt:on # Binarize the already uniformized image else: - LOGGER.info("Stage 2: Skipping skull-strip, generating mask from input") + LOGGER.info("ANAT Stage 2: Skipping skull-strip, generating mask from input") desc += """\ The provided T1w image was previously skull-stripped; a brain mask was derived from the input image. @@ -753,16 +785,21 @@ def init_anat_fit_wf( ]) # fmt:on - ds_mask_wf = init_ds_mask_wf(bids_root=bids_root, output_dir=output_dir) + ds_t1w_mask_wf = init_ds_mask_wf( + bids_root=bids_root, + output_dir=output_dir, + mask_type="brain", + name="ds_t1w_mask_wf", + ) # fmt:off workflow.connect([ - (sourcefile_buffer, ds_mask_wf, [("source_files", "inputnode.source_files")]), - (refined_buffer, ds_mask_wf, [("t1w_mask", "inputnode.t1w_mask")]), - (ds_mask_wf, outputnode, [("outputnode.t1w_mask", "t1w_mask")]), + (sourcefile_buffer, ds_t1w_mask_wf, [("source_files", "inputnode.source_files")]), + (refined_buffer, ds_t1w_mask_wf, [("t1w_mask", "inputnode.mask_file")]), + (ds_t1w_mask_wf, outputnode, [("outputnode.mask_file", "t1w_mask")]), ]) # fmt:on else: - LOGGER.info("Found brain mask") + LOGGER.info("ANAT Found brain mask") desc += """\ A pre-computed brain mask was provided as input and used throughout the workflow. """ @@ -772,7 +809,7 @@ def init_anat_fit_wf( workflow.connect([(anat_validate, apply_mask, [("out_file", "in_file")])]) # Run N4 if it hasn't been pre-run if not have_t1w: - LOGGER.info("Skipping skull-strip, INU-correction only") + LOGGER.info("ANAT Skipping skull-strip, INU-correction only") n4_only_wf = init_n4_only_wf( omp_nthreads=omp_nthreads, atropos_use_random_seed=not skull_strip_fixed_seed, @@ -787,13 +824,13 @@ def init_anat_fit_wf( ]) # fmt:on else: - LOGGER.info("Skipping Stage 2") + LOGGER.info("ANAT Skipping Stage 2") workflow.connect([(apply_mask, t1w_buffer, [("out_file", "t1w_brain")])]) workflow.connect([(refined_buffer, outputnode, [("t1w_mask", "t1w_mask")])]) # Stage 3: Segmentation if not (have_dseg and have_tpms): - LOGGER.info("Stage 3: Preparing segmentation workflow") + LOGGER.info("ANAT Stage 3: Preparing segmentation workflow") fsl_ver = fsl.FAST().version or "(version unknown)" desc += f"""\ Brain tissue segmentation of cerebrospinal fluid (CSF), @@ -833,13 +870,13 @@ def init_anat_fit_wf( ]) # fmt:on else: - LOGGER.info("Skipping Stage 3") + LOGGER.info("ANAT Skipping Stage 3") if have_dseg: - LOGGER.info("Found discrete segmentation") + LOGGER.info("ANAT Found discrete segmentation") desc += "Precomputed discrete tissue segmentations were provided as inputs.\n" seg_buffer.inputs.t1w_dseg = precomputed["t1w_dseg"] if have_tpms: - LOGGER.info("Found tissue probability maps") + LOGGER.info("ANAT Found tissue probability maps") desc += "Precomputed tissue probabiilty maps were provided as inputs.\n" seg_buffer.inputs.t1w_tpms = precomputed["t1w_tpms"] @@ -858,7 +895,7 @@ def init_anat_fit_wf( std2anat_buffer.inputs.in1 = [xfm["reverse"] for xfm in found_xfms.values()] if templates: - LOGGER.info(f"Stage 4: Preparing normalization workflow for {templates}") + LOGGER.info(f"ANAT Stage 4: Preparing normalization workflow for {templates}") register_template_wf = init_register_template_wf( sloppy=sloppy, omp_nthreads=omp_nthreads, @@ -885,7 +922,7 @@ def init_anat_fit_wf( ]) # fmt:on if found_xfms: - LOGGER.info(f"Stage 4: Found pre-computed registrations for {found_xfms}") + LOGGER.info(f"ANAT Stage 4: Found pre-computed registrations for {found_xfms}") # Do not attempt refinement (Stage 6, below) if have_mask or not freesurfer: @@ -901,7 +938,7 @@ def init_anat_fit_wf( workflow.__desc__ = desc if not freesurfer: - LOGGER.info("Skipping Stages 5 and 6") + LOGGER.info("ANAT Skipping Stages 5+") return workflow fs_isrunning = pe.Node( @@ -910,7 +947,7 @@ def init_anat_fit_wf( fs_isrunning.inputs.logger = LOGGER # Stage 5: Surface reconstruction (--fs-no-reconall not set) - LOGGER.info("Stage 5: Preparing surface reconstruction workflow") + LOGGER.info("ANAT Stage 5: Preparing surface reconstruction workflow") surface_recon_wf = init_surface_recon_wf( name="surface_recon_wf", omp_nthreads=omp_nthreads, @@ -956,7 +993,7 @@ def init_anat_fit_wf( ]) # fmt:on elif "reverse" in fsnative_xfms: - LOGGER.info("Found fsnative-T1w transform - skipping registration") + LOGGER.info("ANAT Found fsnative-T1w transform - skipping registration") outputnode.inputs.fsnative2t1w_xfm = fsnative_xfms["reverse"] else: raise RuntimeError( @@ -964,7 +1001,7 @@ def init_anat_fit_wf( ) if not have_mask: - LOGGER.info("Stage 6: Preparing mask refinement workflow") + LOGGER.info("ANAT Stage 6: Preparing mask refinement workflow") # Stage 6: Refine ANTs mask with FreeSurfer segmentation refinement_wf = init_refinement_wf() applyrefined = pe.Node(fsl.ApplyMask(), name="applyrefined") @@ -987,10 +1024,10 @@ def init_anat_fit_wf( ]) # fmt:on else: - LOGGER.info("Found brain mask - skipping Stage 6") + LOGGER.info("ANAT Found brain mask - skipping Stage 6") if t2w and not have_t2w: - LOGGER.info("Stage 7: Creating T2w template") + LOGGER.info("ANAT Stage 7: Creating T2w template") t2w_template_wf = init_anat_template_wf( longitudinal=longitudinal, omp_nthreads=omp_nthreads, @@ -1049,9 +1086,187 @@ def init_anat_fit_wf( ]) # fmt:on elif not t2w: - LOGGER.info("No T2w images provided - skipping Stage 7") + LOGGER.info("ANAT No T2w images provided - skipping Stage 7") + else: + LOGGER.info("ANAT Found preprocessed T2w - skipping Stage 7") + + # Stages 8-10: Surface conversion and registration + # sphere_reg is needed to generate sphere_reg_fsLR + # sphere and sulc are needed to generate sphere_reg_msm + # white, pial, midthickness and thickness are needed to resample in the cortical ribbon + # TODO: Consider paring down or splitting into a subworkflow that can be called on-demand + # A subworkflow would still need to check for precomputed outputs + needed_anat_surfs = ["white", "pial", "midthickness"] + needed_metrics = ["thickness", "sulc"] + needed_spheres = ["sphere_reg", "sphere"] + + # Detect pre-computed surfaces + found_surfs = { + surf: sorted(precomputed[surf]) + for surf in needed_anat_surfs + needed_metrics + needed_spheres + if len(precomputed.get(surf, [])) == 2 + } + if found_surfs: + LOGGER.info(f"ANAT Stage 8: Found pre-converted surfaces for {list(found_surfs)}") + surfaces_buffer.inputs.trait_set(**found_surfs) + + # Stage 8: Surface conversion + surfs = [surf for surf in needed_anat_surfs if surf not in found_surfs] + spheres = [sphere for sphere in needed_spheres if sphere not in found_surfs] + if surfs or spheres: + LOGGER.info(f"ANAT Stage 8: Creating GIFTI surfaces for {surfs + spheres}") + if surfs: + gifti_surfaces_wf = init_gifti_surfaces_wf(surfaces=surfs) + ds_surfaces_wf = init_ds_surfaces_wf( + bids_root=bids_root, output_dir=output_dir, surfaces=surfs + ) + # fmt:off + workflow.connect([ + (surface_recon_wf, gifti_surfaces_wf, [ + ('outputnode.subject_id', 'inputnode.subject_id'), + ('outputnode.subjects_dir', 'inputnode.subjects_dir'), + ('outputnode.fsnative2t1w_xfm', 'inputnode.fsnative2t1w_xfm'), + ]), + (gifti_surfaces_wf, surfaces_buffer, [ + (f'outputnode.{surf}', surf) for surf in surfs + ]), + (sourcefile_buffer, ds_surfaces_wf, [("source_files", "inputnode.source_files")]), + (gifti_surfaces_wf, ds_surfaces_wf, [ + (f'outputnode.{surf}', f'inputnode.{surf}') for surf in surfs + ]), + ]) + # fmt:on + if spheres: + gifti_spheres_wf = init_gifti_surfaces_wf( + surfaces=spheres, to_scanner=False, name='gifti_spheres_wf' + ) + ds_spheres_wf = init_ds_surfaces_wf( + bids_root=bids_root, output_dir=output_dir, surfaces=spheres, name='ds_spheres_wf' + ) + # fmt:off + workflow.connect([ + (surface_recon_wf, gifti_spheres_wf, [ + ('outputnode.subject_id', 'inputnode.subject_id'), + ('outputnode.subjects_dir', 'inputnode.subjects_dir'), + # No transform for spheres, following HCP pipelines' lead + ]), + (gifti_spheres_wf, surfaces_buffer, [ + (f'outputnode.{sphere}', sphere) for sphere in spheres + ]), + (sourcefile_buffer, ds_spheres_wf, [("source_files", "inputnode.source_files")]), + (gifti_spheres_wf, ds_spheres_wf, [ + (f'outputnode.{sphere}', f'inputnode.{sphere}') for sphere in spheres + ]), + ]) + # fmt:on + metrics = [metric for metric in needed_metrics if metric not in found_surfs] + if metrics: + LOGGER.info(f"ANAT Stage 8: Creating GIFTI metrics for {metrics}") + gifti_morph_wf = init_gifti_morphometrics_wf(morphometrics=metrics) + ds_morph_wf = init_ds_surface_metrics_wf( + bids_root=bids_root, output_dir=output_dir, metrics=metrics, name='ds_morph_wf' + ) + + # fmt:off + workflow.connect([ + (surface_recon_wf, gifti_morph_wf, [ + ('outputnode.subject_id', 'inputnode.subject_id'), + ('outputnode.subjects_dir', 'inputnode.subjects_dir'), + ]), + (gifti_morph_wf, surfaces_buffer, [ + (f'outputnode.{metric}', metric) for metric in metrics + ]), + (sourcefile_buffer, ds_morph_wf, [("source_files", "inputnode.source_files")]), + (gifti_morph_wf, ds_morph_wf, [ + (f'outputnode.{metric}', f'inputnode.{metric}') for metric in metrics + ]), + ]) + # fmt:on + + if "anat_ribbon" not in precomputed: + LOGGER.info("ANAT Stage 8a: Creating cortical ribbon mask") + anat_ribbon_wf = init_anat_ribbon_wf() + ds_ribbon_mask_wf = init_ds_mask_wf( + bids_root=bids_root, + output_dir=output_dir, + mask_type="ribbon", + name="ds_ribbon_mask_wf", + ) + # fmt:off + workflow.connect([ + (t1w_buffer, anat_ribbon_wf, [ + ("t1w_preproc", "inputnode.ref_file"), + ]), + (surfaces_buffer, anat_ribbon_wf, [ + ("white", "inputnode.white"), + ("pial", "inputnode.pial"), + ]), + (sourcefile_buffer, ds_ribbon_mask_wf, [("source_files", "inputnode.source_files")]), + (anat_ribbon_wf, ds_ribbon_mask_wf, [ + ("outputnode.anat_ribbon", "inputnode.mask_file"), + ]), + (ds_ribbon_mask_wf, outputnode, [("outputnode.mask_file", "anat_ribbon")]), + ]) + # fmt:on + else: + LOGGER.info("ANAT Stage 8a: Found pre-computed cortical ribbon mask") + outputnode.inputs.anat_ribbon = precomputed["anat_ribbon"] + + # Stage 9: Baseline fsLR registration + if len(precomputed.get("sphere_reg_fsLR", [])) < 2: + LOGGER.info("ANAT Stage 9: Creating fsLR registration sphere") + fsLR_reg_wf = init_fsLR_reg_wf() + ds_fsLR_reg_wf = init_ds_surfaces_wf( + bids_root=bids_root, + output_dir=output_dir, + surfaces=["sphere_reg_fsLR"], + name='ds_fsLR_reg_wf', + ) + + # fmt:off + workflow.connect([ + (surfaces_buffer, fsLR_reg_wf, [('sphere_reg', 'inputnode.sphere_reg')]), + (sourcefile_buffer, ds_fsLR_reg_wf, [("source_files", "inputnode.source_files")]), + (fsLR_reg_wf, ds_fsLR_reg_wf, [ + ('outputnode.sphere_reg_fsLR', 'inputnode.sphere_reg_fsLR') + ]), + (ds_fsLR_reg_wf, fsLR_buffer, [('outputnode.sphere_reg_fsLR', 'sphere_reg_fsLR')]), + ]) + # fmt:on + else: + LOGGER.info("ANAT Stage 9: Found pre-computed fsLR registration sphere") + fsLR_buffer.inputs.sphere_reg_fsLR = sorted(precomputed["sphere_reg_fsLR"]) + + # Stage 10: MSMSulc + if msm_sulc and len(precomputed.get("sphere_reg_msm", [])) < 2: + LOGGER.info("ANAT Stage 10: Creating MSM-Sulc registration sphere") + msm_sulc_wf = init_msm_sulc_wf(sloppy=sloppy) + ds_msmsulc_wf = init_ds_surfaces_wf( + bids_root=bids_root, + output_dir=output_dir, + surfaces=["sphere_reg_msm"], + name='ds_msmsulc_wf', + ) + + # fmt:off + workflow.connect([ + (surfaces_buffer, msm_sulc_wf, [ + ('sulc', 'inputnode.sulc'), + ('sphere', 'inputnode.sphere'), + ]), + (fsLR_buffer, msm_sulc_wf, [('sphere_reg_fsLR', 'inputnode.sphere_reg_fsLR')]), + (sourcefile_buffer, ds_msmsulc_wf, [("source_files", "inputnode.source_files")]), + (msm_sulc_wf, ds_msmsulc_wf, [ + ('outputnode.sphere_reg_fsLR', 'inputnode.sphere_reg_msm') + ]), + (ds_msmsulc_wf, msm_buffer, [('outputnode.sphere_reg_msm', 'sphere_reg_msm')]), + ]) + # fmt:on + elif msm_sulc: + LOGGER.info("ANAT Stage 10: Found pre-computed MSM-Sulc registration sphere") + fsLR_buffer.inputs.sphere_reg_msm = sorted(precomputed["sphere_reg_msm"]) else: - LOGGER.info("Found preprocessed T2w - skipping Stage 7") + LOGGER.info("ANAT Stage 10: MSM-Sulc disabled") return workflow diff --git a/smriprep/workflows/base.py b/smriprep/workflows/base.py index 8769673c60..008c172e58 100644 --- a/smriprep/workflows/base.py +++ b/smriprep/workflows/base.py @@ -366,7 +366,7 @@ def init_single_subject_wf( std_spaces = spaces.get_spaces(nonstandard=False, dim=(3,)) std_spaces.append("fsnative") for deriv_dir in derivatives: - deriv_cache.update(collect_derivatives(deriv_dir, subject_id, std_spaces, freesurfer)) + deriv_cache.update(collect_derivatives(deriv_dir, subject_id, std_spaces)) inputnode = pe.Node(niu.IdentityInterface(fields=["subjects_dir"]), name="inputnode") diff --git a/smriprep/workflows/outputs.py b/smriprep/workflows/outputs.py index ecefed25f5..99211b0faa 100644 --- a/smriprep/workflows/outputs.py +++ b/smriprep/workflows/outputs.py @@ -316,8 +316,9 @@ def init_ds_template_wf( def init_ds_mask_wf( *, - bids_root, - output_dir, + bids_root: str, + output_dir: str, + mask_type: str, name="ds_mask_wf", ): """ @@ -336,40 +337,48 @@ def init_ds_mask_wf( ------ source_files List of input T1w images - t1w_mask - Mask of the ``t1w_preproc`` + mask_file + Mask to save Outputs ------- - t1w_mask - The location in the output directory of the T1w mask + mask_file + The location in the output directory of the mask """ workflow = Workflow(name=name) inputnode = pe.Node( - niu.IdentityInterface(fields=["source_files", "t1w_mask"]), + niu.IdentityInterface(fields=["source_files", "mask_file"]), name="inputnode", ) - outputnode = pe.Node(niu.IdentityInterface(fields=["t1w_mask"]), name="outputnode") + outputnode = pe.Node(niu.IdentityInterface(fields=["mask_file"]), name="outputnode") raw_sources = pe.Node(niu.Function(function=_bids_relative), name="raw_sources") raw_sources.inputs.bids_root = bids_root - ds_t1w_mask = pe.Node( - DerivativesDataSink(base_directory=output_dir, desc="brain", suffix="mask", compress=True), + ds_mask = pe.Node( + DerivativesDataSink( + base_directory=output_dir, + desc=mask_type, + suffix="mask", + compress=True, + ), name="ds_t1w_mask", run_without_submitting=True, ) - ds_t1w_mask.inputs.Type = "Brain" + if mask_type == "brain": + ds_mask.inputs.Type = "Brain" + else: + ds_mask.inputs.Type = "ROI" # fmt:off workflow.connect([ (inputnode, raw_sources, [('source_files', 'in_files')]), - (inputnode, ds_t1w_mask, [('t1w_mask', 'in_file'), - ('source_files', 'source_file')]), - (raw_sources, ds_t1w_mask, [('out', 'RawSources')]), - (ds_t1w_mask, outputnode, [('out_file', 't1w_mask')]), + (inputnode, ds_mask, [('mask_file', 'in_file'), + ('source_files', 'source_file')]), + (raw_sources, ds_mask, [('out', 'RawSources')]), + (ds_mask, outputnode, [('out_file', 'mask_file')]), ]) # fmt:on @@ -722,10 +731,7 @@ def init_ds_surfaces_wf( # fmt:off workflow.connect([ - (inputnode, ds_surf, [ - (surf, 'in_file'), - ('source_files', 'source_file'), - ]), + (inputnode, ds_surf, [(surf, 'in_file'), ('source_files', 'source_file')]), (ds_surf, outputnode, [('out_file', surf)]), ]) # fmt:on @@ -733,6 +739,71 @@ def init_ds_surfaces_wf( return workflow +def init_ds_surface_metrics_wf( + *, + bids_root: str, + output_dir: str, + metrics: list[str], + name="ds_surface_metrics_wf", +) -> Workflow: + """ + Save GIFTI surface metrics + + Parameters + ---------- + bids_root : :class:`str` + Root path of BIDS dataset + output_dir : :class:`str` + Directory in which to save derivatives + metrics : :class:`str` + List of metrics to generate DataSinks for + name : :class:`str` + Workflow name (default: ds_surface_metrics_wf) + + Inputs + ------ + source_files + List of input T1w images + ```` + Left and right GIFTIs for each metric passed to ``metrics`` + + Outputs + ------- + ```` + Left and right GIFTIs in ``output_dir`` for each metric passed to ``metrics`` + + """ + workflow = Workflow(name=name) + + inputnode = pe.Node( + niu.IdentityInterface(fields=["source_files"] + metrics), + name="inputnode", + ) + outputnode = pe.Node(niu.IdentityInterface(fields=metrics), name="outputnode") + + for metric in metrics: + ds_surf = pe.MapNode( + DerivativesDataSink( + base_directory=output_dir, + hemi=["L", "R"], + suffix=metric, + extension=".shape.gii", + ), + iterfield=("in_file", "hemi"), + name=f"ds_{metric}", + run_without_submitting=True, + ) + + # fmt:off + workflow.connect([ + (inputnode, ds_surf, [(metric, 'in_file'), ('source_files', 'source_file')]), + (ds_surf, outputnode, [('out_file', metric)]), + ]) + # fmt:on + + return workflow + + def init_anat_second_derivatives_wf( *, bids_root: str, @@ -803,11 +874,8 @@ def init_anat_second_derivatives_wf( "t1w_dseg", "t1w_tpms", "anat2std_xfm", - "surfaces", "sphere_reg", "sphere_reg_fsLR", - "morphometrics", - "anat_ribbon", "t1w_fs_aseg", "t1w_fs_aparc", "cifti_morph", @@ -942,34 +1010,6 @@ def init_anat_second_derivatives_wf( if not freesurfer: return workflow - from niworkflows.interfaces.surf import Path2BIDS - - # Morphometrics - name_morphs = pe.MapNode( - Path2BIDS(), - iterfield="in_file", - name="name_morphs", - run_without_submitting=True, - ) - ds_morphs = pe.MapNode( - DerivativesDataSink(base_directory=output_dir, extension=".shape.gii"), - iterfield=["in_file", "hemi", "suffix"], - name="ds_morphs", - run_without_submitting=True, - ) - # Ribbon volume - ds_anat_ribbon = pe.Node( - DerivativesDataSink( - base_directory=output_dir, - desc="ribbon", - suffix="mask", - extension=".nii.gz", - compress=True, - ), - name="ds_anat_ribbon", - run_without_submitting=True, - ) - # Parcellations ds_t1w_fsaseg = pe.Node( DerivativesDataSink(base_directory=output_dir, desc="aseg", suffix="dseg", compress=True), @@ -986,18 +1026,10 @@ def init_anat_second_derivatives_wf( # fmt:off workflow.connect([ - (inputnode, name_morphs, [('morphometrics', 'in_file')]), - (inputnode, ds_morphs, [('morphometrics', 'in_file'), - ('source_files', 'source_file')]), - (name_morphs, ds_morphs, [('hemi', 'hemi'), - ('suffix', 'suffix')]), (inputnode, ds_t1w_fsaseg, [('t1w_fs_aseg', 'in_file'), ('source_files', 'source_file')]), (inputnode, ds_t1w_fsparc, [('t1w_fs_aparc', 'in_file'), ('source_files', 'source_file')]), - (inputnode, ds_anat_ribbon, [('anat_ribbon', 'in_file'), - ('source_files', 'source_file')]), - ]) # fmt:on diff --git a/smriprep/workflows/surfaces.py b/smriprep/workflows/surfaces.py index b682b53cb1..90a885a0c7 100644 --- a/smriprep/workflows/surfaces.py +++ b/smriprep/workflows/surfaces.py @@ -579,9 +579,7 @@ def _dedup(in_list): def init_surface_derivatives_wf( *, - msm_sulc: bool = False, cifti_output: ty.Literal["91k", "170k", False] = False, - sloppy: bool = False, name="surface_derivatives_wf", ): r""" @@ -639,15 +637,8 @@ def init_surface_derivatives_wf( outputnode = pe.Node( niu.IdentityInterface( fields=[ - "surfaces", - "white", - "pial", - "midthickness", "inflated", - "morphometrics", - "sphere_reg", - "sphere_reg_fsLR", - "sphere_reg_msm", + "curv", "out_aseg", "out_aparc", "cifti_morph", @@ -657,15 +648,8 @@ def init_surface_derivatives_wf( name="outputnode", ) - gifti_surfaces_wf = init_gifti_surfaces_wf() - gifti_spheres_wf = init_gifti_surfaces_wf( - surfaces=["sphere", "sphere_reg"], - to_scanner=False, - name="gifti_spheres_wf", - ) - gifti_morph_wf = init_gifti_morphometrics_wf() - fsLR_reg_wf = init_fsLR_reg_wf() - msm_sulc_wf = init_msm_sulc_wf(sloppy=sloppy) + gifti_surfaces_wf = init_gifti_surfaces_wf(surfaces=["inflated"]) + gifti_morph_wf = init_gifti_morphometrics_wf(morphometrics=["curv"]) aseg_to_native_wf = init_segs_to_native_wf() aparc_to_native_wf = init_segs_to_native_wf(segmentation="aparc_aseg") @@ -677,22 +661,10 @@ def init_surface_derivatives_wf( ('subject_id', 'inputnode.subject_id'), ('fsnative2t1w_xfm', 'inputnode.fsnative2t1w_xfm'), ]), - (inputnode, gifti_spheres_wf, [ - ('subjects_dir', 'inputnode.subjects_dir'), - ('subject_id', 'inputnode.subject_id'), - ]), (inputnode, gifti_morph_wf, [ ('subjects_dir', 'inputnode.subjects_dir'), ('subject_id', 'inputnode.subject_id'), ]), - (gifti_spheres_wf, fsLR_reg_wf, [ - ('outputnode.sphere_reg', 'inputnode.sphere_reg'), - ]), - (gifti_spheres_wf, msm_sulc_wf, [('outputnode.sphere', 'inputnode.sphere')]), - (fsLR_reg_wf, msm_sulc_wf, [ - ('outputnode.sphere_reg_fsLR', 'inputnode.sphere_reg_fsLR'), - ]), - (gifti_morph_wf, msm_sulc_wf, [('outputnode.sulc', 'inputnode.sulc')]), (inputnode, aseg_to_native_wf, [ ('subjects_dir', 'inputnode.subjects_dir'), ('subject_id', 'inputnode.subject_id'), @@ -707,19 +679,10 @@ def init_surface_derivatives_wf( ]), # Output - (gifti_surfaces_wf, outputnode, [ - ('outputnode.surfaces', 'surfaces'), - ('outputnode.white', 'white'), - ('outputnode.pial', 'pial'), - ('outputnode.midthickness', 'midthickness'), - ('outputnode.inflated', 'inflated'), - ]), - (gifti_spheres_wf, outputnode, [('outputnode.sphere_reg', 'sphere_reg')]), - (gifti_morph_wf, outputnode, [('outputnode.morphometrics', 'morphometrics')]), - (fsLR_reg_wf, outputnode, [('outputnode.sphere_reg_fsLR', 'sphere_reg_fsLR')]), - (msm_sulc_wf, outputnode, [('outputnode.sphere_reg_fsLR', 'sphere_reg_msm')]), + (gifti_surfaces_wf, outputnode, [('outputnode.inflated', 'inflated')]), (aseg_to_native_wf, outputnode, [('outputnode.out_file', 'out_aseg')]), (aparc_to_native_wf, outputnode, [('outputnode.out_file', 'out_aparc')]), + (gifti_morph_wf, outputnode, [('outputnode.curv', 'curv')]), ]) # fmt:on @@ -1004,7 +967,7 @@ def init_gifti_morphometrics_wf( workflow = Workflow(name=name) inputnode = pe.Node( - niu.IdentityInterface(["subjects_dir", "subject_id", "fsnative2t1w_xfm"]), + niu.IdentityInterface(["subjects_dir", "subject_id"]), name="inputnode", ) outputnode = pe.Node(