From 5accc40ba2defc0cac13510b09526b12608e6455 Mon Sep 17 00:00:00 2001 From: Ali Khan Date: Sun, 29 Dec 2024 13:00:33 -0500 Subject: [PATCH] close to final working version - only one atlas from CLI instead of multiples (was using atlas[0] some places) - remove final transforms (not used with nativesurf workflow) - giftis only for atlas (needs updated multihist7 files - workflow to generate these in resources/atlas-v2 for now) - use unfoldedreg when the atlas has the metrics - option for dentate-based unfoldreg too (unused right now) TODO: - update atlas files (use another location for testing) - check unfoldreg warps - make nativesurf the default --- hippunfold/config/snakebids.yml | 26 +- hippunfold/resources/atlas-v2/README.md | 7 + hippunfold/resources/atlas-v2/Snakefile | 90 +++++ .../resources/atlas-v2/labellist_withdg.txt | 12 + hippunfold/workflow/Snakefile | 4 - hippunfold/workflow/rules/common.smk | 70 +--- hippunfold/workflow/rules/gifti.smk | 95 +---- hippunfold/workflow/rules/native_surf.smk | 335 ++++++++++++------ .../rules/resample_final_to_crop_native.smk | 2 +- hippunfold/workflow/rules/subfields.smk | 78 +--- hippunfold/workflow/rules/unfoldreg_v1.smk | 48 +++ hippunfold/workflow/rules/warps.smk | 230 ------------ 12 files changed, 422 insertions(+), 575 deletions(-) create mode 100644 hippunfold/resources/atlas-v2/README.md create mode 100755 hippunfold/resources/atlas-v2/Snakefile create mode 100755 hippunfold/resources/atlas-v2/labellist_withdg.txt diff --git a/hippunfold/config/snakebids.yml b/hippunfold/config/snakebids.yml index 850a729d..abb08af5 100644 --- a/hippunfold/config/snakebids.yml +++ b/hippunfold/config/snakebids.yml @@ -187,11 +187,10 @@ parse_args: - 'mouse' default: - 'multihist7' - nargs: '+' help: 'Select the atlas (unfolded space) to use for subfield labels. (default: %(default)s)' --no_unfolded_reg: - help: 'Do not perform unfolded space (2D) registration based on thickness, curvature, and gyrification for closer alignment to the reference atlas. NOTE: only multihist7 has these features currently, so this unfolded_reg is automatically skipped if a different atlas is chosen. (default: %(default)s)' + help: 'Do not perform unfolded space (2D) registration based on surface metrics (e.g. thickness, curvature, and gyrification) for closer alignment to the reference atlas. NOTE: only multihist7 has these features currently, so this unfolded_reg is automatically skipped if a different atlas is chosen. (default: %(default)s)' default: False action: store_true @@ -450,12 +449,21 @@ template_files: atlas_files: multihist7: - label_nii: sub-maxprob_label-hipp_desc-manualsubfieldsunfoldaligned_dseg.nii.gz - lut: desc-subfields_atlas-multihist7_dseg.tsv - label_list: labellist.txt - thickness: thickness.nii.gz - curvature: curvature.nii.gz - gyrification: gyrification.nii.gz + lut: desc-subfields_atlas-multihist7_dseg.tsv #<-- perhaps generate this from the cifti instead of having user provide it + surf_gii: hemi-{hemi}_space-unfold_label-{label}_{surf_type}.surf.gii + metric_gii: hemi-{hemi}_label-{label}_{metric}.shape.gii + label_gii: hemi-{hemi}_label-{label}_dseg.label.gii + label_wildcards: + - hipp + hemi_wildcards: + - L + - R + metric_wildcards: + - thickness + - gyrification + - curvature + + bigbrain: label_nii: sub-bigbrain_hemi-{hemi}_label-hipp_desc-manualsubfields_dseg.nii.gz lut: desc-subfields_atlas-bigbrain_dseg.tsv @@ -536,6 +544,7 @@ resource_urls: synthseg_v0.1: 'zenodo.org/record/8184230/files/trained_model.3d_fullres.Task102_synsegGenDetailed.nnUNetTrainerV2.model_best.tar' synthseg_v0.2: 'zenodo.org/record/8184230/files/trained_model.3d_fullres.Task203_synthseg.nnUNetTrainerV2.model_best.tar' atlas: + fakeatlas: 'this.is.a/fake/url' multihist7: 'files.ca-1.osf.io/v1/resources/v8acf/providers/osfstorage/65395b782827451220b86dd8/?zip=' bigbrain: 'files.ca-1.osf.io/v1/resources/v8acf/providers/osfstorage/65395b8b13d27b123094c96f/?zip=' magdeburg: 'files.ca-1.osf.io/v1/resources/v8acf/providers/osfstorage/65395b8013d27b122f94c938/?zip=' @@ -691,7 +700,6 @@ skip_inject_template_labels: False force_nnunet_model: False t1_reg_template: False generate_myelin_map: False -no_unfolded_reg: False root: results use_template_seg: False template_seg_smoothing_factor: 2 diff --git a/hippunfold/resources/atlas-v2/README.md b/hippunfold/resources/atlas-v2/README.md new file mode 100644 index 00000000..4fc579a9 --- /dev/null +++ b/hippunfold/resources/atlas-v2/README.md @@ -0,0 +1,7 @@ +# WIP - workflow for updating atlases to new v2 specs + +This example is run from the existing multihist atlas folder. + +Still needs updates for cifti structure labels etc (will add this in once we have an updated wb_command/wb_view in our deps) + +We currently have templates and atlases in hippunfold (selected by --template and --atlas), where the former selects the volume-based template used for volumetric registration, and the latter selects the unfolded-space subfields (and metrics for unfolded registration). I am proposing the latter become purely surface-based, so instead of providing labels and metrics in the volumetric unfolded space, we use giftis for everything, and this now becomes surface-templates instead of atlases. diff --git a/hippunfold/resources/atlas-v2/Snakefile b/hippunfold/resources/atlas-v2/Snakefile new file mode 100755 index 00000000..b4ee4473 --- /dev/null +++ b/hippunfold/resources/atlas-v2/Snakefile @@ -0,0 +1,90 @@ + +hippunfold_dir='/localscratch/hippunfold_khanlab' +hipp_dseg='sub-maxprob_label-hipp_desc-manualsubfieldsunfoldaligned_dseg.nii.gz' +labellist='labellist_withdg.txt' + +hemis=['L','R'] +metrics=['thickness','curvature','gyrification'] +surftypes=['midthickness','inner','outer'] + +hemi_to_structure = {"L": "CORTEX_LEFT", "R": "CORTEX_RIGHT"} +surf_to_secondary_type = { + "midthickness": "MIDTHICKNESS", + "inner": "PIAL", + "outer": "GRAY_WHITE", +} + + + +rule all: + input: + expand('hemi-{hemi}_label-{label}_dseg.label.gii',label=['hipp','dentate'],hemi=hemis), + expand('hemi-{hemi}_label-{label}_{metric}.shape.gii',label=['hipp'],hemi=hemis,metric=metrics), + expand('hemi-{hemi}_space-unfold_label-{label}_{surftype}.surf.gii',label=['hipp','dentate'],hemi=hemis,surftype=surftypes), + + + +rule hipp_dseg: + input: + hipp_dseg + output: + temp('hipp_dseg.nii.gz') + shell: + 'cp {input} {output}' + +rule dentate_dseg: + input: + hipp_dseg + params: + dg_label=6 + output: + temp('dentate_dseg.nii.gz') + shell: + 'c3d {input} -scale 0 -shift {params.dg_label} -o {output}' + +rule surf_gifti: + input: + surf_gii=Path(hippunfold_dir) / 'hippunfold' / 'resources' / 'unfold_template_{label}' / 'tpl-avg_space-unfold_den-unfoldiso_{surftype}.surf.gii' + params: + structure_type=lambda wildcards: hemi_to_structure[wildcards.hemi], + secondary_type=lambda wildcards: surf_to_secondary_type[wildcards.surftype], + surface_type="FLAT", + output: + surf_gii='hemi-{hemi}_space-unfold_label-{label}_{surftype}.surf.gii' + shell: + 'cp {input} {output} && ' + "wb_command -set-structure {output} {params.structure_type} -surface-type {params.surface_type}" + " -surface-secondary-type {params.secondary_type}" + + +rule label_gifti: + input: + dseg_nii='{label}_dseg.nii.gz', + lut=labellist, + surf_gii='hemi-{hemi}_space-unfold_label-{label}_midthickness.surf.gii' + params: + structure_type=lambda wildcards: hemi_to_structure[wildcards.hemi], + output: + label_gii='hemi-{hemi}_label-{label}_dseg.label.gii' + shadow: 'minimal' + shell: + 'wb_command -volume-label-import {input.dseg_nii} {input.lut} label_vol.nii && ' + 'wb_command -volume-label-to-surface-mapping label_vol.nii {input.surf_gii} {output.label_gii} && ' + "wb_command -set-structure {output} {params.structure_type}" + + + +rule metric_gifti: + input: + metric_nii='{metric}.nii.gz', + surf_gii='hemi-{hemi}_space-unfold_label-{label}_inner.surf.gii' #existing metric niftis are on inner space + params: + structure_type=lambda wildcards: hemi_to_structure[wildcards.hemi], + output: + metric_gii='hemi-{hemi}_label-{label}_{metric}.shape.gii' + shadow: 'minimal' + shell: + 'wb_command -volume-to-surface-mapping {input.metric_nii} {input.surf_gii} {output.metric_gii} -trilinear && ' + "wb_command -set-structure {output} {params.structure_type}" + + diff --git a/hippunfold/resources/atlas-v2/labellist_withdg.txt b/hippunfold/resources/atlas-v2/labellist_withdg.txt new file mode 100755 index 00000000..1aecd491 --- /dev/null +++ b/hippunfold/resources/atlas-v2/labellist_withdg.txt @@ -0,0 +1,12 @@ +Subiculum +1 0 0 255 255 +CA1 +2 133 222 255 255 +CA2 +3 0 255 170 255 +CA3 +4 255 162 0 255 +CA4 +5 255 0 0 255 +DG +6 255 255 0 255 diff --git a/hippunfold/workflow/Snakefile b/hippunfold/workflow/Snakefile index 806ff004..d4f5389a 100644 --- a/hippunfold/workflow/Snakefile +++ b/hippunfold/workflow/Snakefile @@ -61,10 +61,6 @@ if config["skip_inject_template_labels"]: # TODO do we need this? limit_to_list.add(config["modality"]) -# if atlas doesn't contain alignment features -if not config["atlas"] == ["multihist7"]: - config["no_unfolded_reg"] = True - inputs = snakebids.generate_inputs( bids_dir=config["bids_dir"], diff --git a/hippunfold/workflow/rules/common.smk b/hippunfold/workflow/rules/common.smk index 64a20e1e..a6e78db8 100644 --- a/hippunfold/workflow/rules/common.smk +++ b/hippunfold/workflow/rules/common.smk @@ -155,67 +155,6 @@ def get_final_coords(): return coords -def get_final_transforms(): - xfms = [] - - xfms.extend( - inputs[config["modality"]].expand( - bids( - root=root, - datatype="warps", - **inputs.subj_wildcards, - label="{autotop}", - suffix="xfm.nii.gz", - hemi="{hemi}", - from_="{space}", - to="unfold", - mode="image", - ), - space=ref_spaces, - autotop=config["autotop_labels"], - hemi=config["hemi"], - allow_missing=True, - ) - ) - - xfms.extend( - inputs[config["modality"]].expand( - bids( - root=root, - datatype="warps", - **inputs.subj_wildcards, - label="{autotop}", - suffix="xfm.nii.gz", - hemi="{hemi}", - from_="unfold", - to="{space}", - mode="image", - ), - space=ref_spaces, - autotop=config["autotop_labels"], - hemi=config["hemi"], - allow_missing=True, - ) - ) - - xfms.extend( - inputs[config["modality"]].expand( - bids( - root=root, - datatype="warps", - **inputs.subj_wildcards, - label="{autotop}", - suffix="refvol.nii.gz", - space="unfold", - ), - autotop=config["autotop_labels"], - allow_missing=True, - ) - ) - - return xfms - - def get_final_anat(): anat = [] if "T1w" in ref_spaces or "T2w" in ref_spaces: @@ -335,11 +274,10 @@ def get_final_qc(): def get_final_subj_output(): subj_output = [] subj_output.extend(get_final_spec()) - # subj_output.extend(get_final_subfields()) - # subj_output.extend(get_final_coords()) - # subj_output.extend(get_final_transforms()) - # subj_output.extend(get_final_anat()) - # subj_output.extend(get_final_qc()) + subj_output.extend(get_final_subfields()) + subj_output.extend(get_final_coords()) + subj_output.extend(get_final_anat()) + subj_output.extend(get_final_qc()) return subj_output diff --git a/hippunfold/workflow/rules/gifti.smk b/hippunfold/workflow/rules/gifti.smk index 5fd3cc34..5384a4f0 100644 --- a/hippunfold/workflow/rules/gifti.smk +++ b/hippunfold/workflow/rules/gifti.smk @@ -214,98 +214,6 @@ rule affine_gii_to_native: shell: "wb_command -surface-apply-affine {input.gii} {input.xfm} {output.gii}" -rule resample_atlas_to_refvol: - """this is just done in case the atlas has a different unfolded config than the current run""" - input: - refvol=bids( - root=root, - space="unfold", - label="hipp", - datatype="warps", - suffix="refvol.nii.gz", - **inputs.subj_wildcards - ), - atlas_dir=lambda wildcards: Path(download_dir) / "atlas" / wildcards.atlas, - params: - atlas=lambda wildcards, input: Path(input.atlas_dir) - / config["atlas_files"][wildcards.atlas]["label_nii"].format(**wildcards), - output: - label_nii=bids( - root=work, - datatype="anat", - suffix="subfields.nii.gz", - space="unfold", - hemi="{hemi}", - label="hipp", - atlas="{atlas}", - **inputs.subj_wildcards - ), - log: - bids( - root="logs", - suffix="resamplesubfieldrefvol", - space="unfold", - hemi="{hemi}", - label="hipp", - atlas="{atlas}", - **inputs.subj_wildcards - ), - container: - config["singularity"]["autotop"] - group: - "subj" - shell: - "antsApplyTransforms -d 3 -n MultiLabel -i {params.atlas} -r {input.refvol} -o {output.label_nii} -v &> {log}" - - -rule nii_to_label_gii: - input: - label_nii=bids( - root=work, - datatype="anat", - suffix="subfields.nii.gz", - space="unfold", - hemi="{hemi}", - label="hipp", - atlas="{atlas}", - **inputs.subj_wildcards - ), - surf=os.path.join( - workflow.basedir, - "..", - "resources", - "unfold_template_hipp", - "tpl-avg_space-unfold_den-{density}_midthickness.surf.gii", - ), - atlas_dir=lambda wildcards: Path(download_dir) / "atlas" / wildcards.atlas, - params: - label_list=lambda wildcards, input: Path(input.atlas_dir) - / config["atlas_files"][wildcards.atlas]["label_list"].format(**wildcards), - structure_type=lambda wildcards: hemi_to_structure[wildcards.hemi], - output: - label_gii=bids( - root=root, - datatype="surf", - den="{density}", - suffix="subfields.label.gii", - space="{space}", - hemi="{hemi}", - label="hipp", - atlas="{atlas}", - **inputs.subj_wildcards - ), - group: - "subj" - container: - config["singularity"]["autotop"] - shadow: - "minimal" - shell: - "wb_command -volume-to-surface-mapping {input.label_nii} {input.surf} temp.shape.gii -enclosing && " - "wb_command -metric-label-import temp.shape.gii {params.label_list} {output.label_gii} && " - "wb_command -set-structure {output.label_gii} {params.structure_type}" - - def get_cmd_cifti_metric(wildcards, input, output): cmd = f"wb_command -cifti-create-dense-scalar {output}" if "L" in config["hemi"]: @@ -505,8 +413,7 @@ rule create_spec_file_hipp: **inputs.subj_wildcards ), surfname=["midthickness"], - space="{space}", - #space=["{space}", "unfold"], + space=["{space}", "unfold","unfoldreg"], allow_missing=True, ), cifti_metrics=lambda wildcards: inputs[config["modality"]].expand( diff --git a/hippunfold/workflow/rules/native_surf.smk b/hippunfold/workflow/rules/native_surf.smk index 81a6b418..f456d7f5 100644 --- a/hippunfold/workflow/rules/native_surf.smk +++ b/hippunfold/workflow/rules/native_surf.smk @@ -1,3 +1,13 @@ +# This is one big rules file for now, but could split it up according to the +# dividers ( --- ) that describe each conceptual part of the pipeline + +# --- parameters - will put these in config eventually + +wildcard_constraints: + label='[0-9a-zA-Z]+', + metric='[0-9a-zA-Z]+', + native_modality='[0-9a-zA-Z]+' + surf_thresholds = {"inner": 0, "outer": 1, "midthickness": 0.5} # this is for the mapping from inner to outer @@ -20,11 +30,12 @@ resample_from_density = ( unfoldreg_metrics = [ 'thickness', 'curvature', 'gyrification'] -# if we aren't doing unfolded reg, then we have these rules override correct_bad_vertices2 -ruleorder: resample_native_surf > cp_template_to_unfold -#ruleorder: resample_native_surf > correct_bad_vertices2 -#ruleorder: warp_midthickness_to_inout > correct_bad_vertices2 +ruleorder: resample_native_surf_to_std_density > cp_template_to_unfold ruleorder: unfold_spring_model > warp_unfold_native_to_unfoldreg +ruleorder: atlas_label_to_unfold_nii > atlas_metric_to_unfold_nii #temporary until we change from subfields.nii to desc-subfields_dseg.nii for unfold space + +#--- isosurface generation --- + rule prep_hipp_coords_for_meshing: input: @@ -229,6 +240,7 @@ rule update_native_mesh_structure: "cp {input} {output} && wb_command -set-structure {output.surf_gii} {params.structure_type} -surface-type {params.surface_type}" " -surface-secondary-type {params.secondary_type}" +#--- creating unfold surface from native anatomical rule warp_native_mesh_to_unfold: input: @@ -275,6 +287,9 @@ rule warp_native_mesh_to_unfold: "wb_command -set-structure {output.surf_gii} {params.structure_type} -surface-type {params.surface_type}" " -surface-secondary-type {params.secondary_type}" + +# --- creating inner/outer surfaces from native anatomical (using 3d label deformable registration) + rule compute_halfthick_mask: input: coords=bids( @@ -497,6 +512,8 @@ rule warp_midthickness_to_inout: " -surface-secondary-type {params.secondary_type}" +# --- affine transforming anatomical surfaces from corobl to other (T1w, T2w) spaces + # warp native surface from corobl to T1w/T2w rule affine_gii_corobl_to_modality: input: @@ -538,6 +555,8 @@ rule affine_gii_corobl_to_modality: +# --- WIP rule for adjusting mesh in unfolded space + rule unfold_spring_model: input: surf_gii=bids( @@ -566,6 +585,8 @@ rule unfold_spring_model: script: '../scripts/unfold_spring_model.py' +# --- calculating metrics on native anatomical surfaces + rule calculate_surface_area: input: @@ -704,6 +725,7 @@ rule calculate_thickness_from_surface: "wb_command -surface-to-surface-3d-distance {input.outer} {input.inner} {output}" +# --- unfolded registration #for unfold_reg, need to get data into 2d image - Jordan's previous code used metric-to-volume-mapping with unfoldiso inner and outer surfaces @@ -716,10 +738,14 @@ rule pad_unfold_ref: to get a 2D slice we will take the central slice (the next rule creates the single-slice reference for this).""" input: - atlas_dir=Path(download_dir) / "atlas" / "multihist7", - params: - ref_nii=lambda wildcards, input: Path(input.atlas_dir) - / config["atlas_files"]["multihist7"]["label_nii"].format(**wildcards), + ref_nii=bids( + root=root, + space="unfold", + label="{label}", + datatype="warps", + suffix="refvol.nii.gz", + **inputs.subj_wildcards + ), output: ref_nii=bids( root=work, @@ -736,8 +762,8 @@ rule pad_unfold_ref: group: "subj" shell: - #"c3d {params.ref_nii} -scale 0 -shift 1 -pad 16x16x0vox 16x16x0vox 0 -o {output.ref_nii}" - "c3d {params.ref_nii} -scale 0 -shift 1 -pad 64x64x0vox 64x64x0vox 0 -o {output.ref_nii}" + #"c3d {input.ref_nii} -scale 0 -shift 1 -pad 16x16x0vox 16x16x0vox 0 -o {output.ref_nii}" + "c3d {input.ref_nii} -scale 0 -shift 1 -pad 64x64x0vox 64x64x0vox 0 -o {output.ref_nii}" rule extract_unfold_ref_slice: @@ -771,6 +797,8 @@ rule extract_unfold_ref_slice: shell: "c3d {input.ref_3d_nii} -slice z 50% -o {output.ref_2d_nii}" + +#TODO - if we add density wildcard here then it will make use of the template unfold standard densities rule native_metric_to_unfold_nii: """converts metric .gii files to .nii for use in ANTs""" input: @@ -841,12 +869,56 @@ rule native_metric_to_unfold_nii: " -ribbon-constrained {input.inner_surf} {input.outer_surf}" - -rule unfoldreg: - """performs 2D registration in unfolded space as in [reference paper] - this is done in a shadow directory to get rid of the tmp files generated by ants.""" + +#TODO - if we add density wildcard for inputs here, then it will make use of the template unfold standard densities, +# if no density wildcard, then will use the native mesh +# +# native metric, resampled to a standard density + +rule atlas_metric_to_unfold_nii: + """converts metric .gii files to .nii for use in ANTs. + This rule is for the surface template""" input: - fixed_images=expand(bids( + atlas_dir=lambda wildcards: Path(download_dir) / "atlas" / wildcards.atlas, + ref_nii=bids( + root=work, + datatype="anat", + suffix="refvol.nii.gz", + space="unfold", + desc="slice", + hemi="{hemi}", + label="{label}", + **inputs.subj_wildcards + ), #we use ref_nii defined in this workflow so that parameters like spacing and padding of unfolded space can be defined on the fly + params: + metric_gii = lambda wildcards, input: Path(input.atlas_dir) / config["atlas_files"][wildcards.atlas]['metric_gii'].format(**wildcards), + inner_surf=lambda wildcards, input: Path(input.atlas_dir) / config["atlas_files"][wildcards.atlas]['surf_gii'].format(surf_type='inner',**wildcards), + outer_surf=lambda wildcards, input: Path(input.atlas_dir) / config["atlas_files"][wildcards.atlas]['surf_gii'].format(surf_type='outer',**wildcards), + midthickness_surf=lambda wildcards, input: Path(input.atlas_dir) / config["atlas_files"][wildcards.atlas]['surf_gii'].format(surf_type='midthickness',**wildcards), + output: + metric_nii=bids( + root=work, + datatype="anat", + suffix="{metric}.nii.gz", + space="unfold", + hemi="{hemi}", + label="{label}", + atlas='{atlas}', + **inputs.subj_wildcards + ), + container: + config["singularity"]["autotop"] + group: + "subj" + shell: + "wb_command -metric-to-volume-mapping {params.metric_gii} {params.midthickness_surf} {input.ref_nii} {output.metric_nii} " + " -ribbon-constrained {params.inner_surf} {params.outer_surf}" + + +def get_fixed_images_unfoldreg(wildcards): + unfoldreg_metrics = config['atlas_files'][wildcards.atlas]['metric_wildcards'] + + return expand(bids( root=work, datatype="anat", suffix="{metric}.nii.gz", @@ -855,15 +927,36 @@ rule unfoldreg: label="{label}", **inputs.subj_wildcards ), - metric=unfoldreg_metrics,allow_missing=True), - atlas_dir=lambda wildcards: Path(download_dir) / "atlas" / wildcards.atlas, + metric=unfoldreg_metrics,**wildcards) + +def get_moving_images_unfoldreg(wildcards): + unfoldreg_metrics = config['atlas_files'][wildcards.atlas]['metric_wildcards'] + return expand(bids( + root=work, + datatype="anat", + suffix="{metric}.nii.gz", + space="unfold", + hemi="{hemi}", + label="{label}", + atlas="{atlas}", + **inputs.subj_wildcards + ), + metric=unfoldreg_metrics,**wildcards) + + + +rule unfoldreg: + """performs 2D registration in unfolded space as in [reference paper] + this is done in a shadow directory to get rid of the tmp files generated by ants.""" + input: + fixed_images=get_fixed_images_unfoldreg, + moving_images=get_moving_images_unfoldreg params: antsparams="-d 2 -t so -o tmp", fixed_args=lambda wildcards, input: ' '.join(['-f {img}'.format(img=img) for img in input.fixed_images]), moving_args=lambda wildcards, input: - ' '.join(['-m {img}'.format(img=str( Path(input.atlas_dir) / config["atlas_files"][wildcards.atlas][metric].format(**wildcards))) - for metric in unfoldreg_metrics ] ), + ' '.join(['-m {img}'.format(img=img) for img in input.moving_images]), cmd_copy_warps=lambda wildcards, output: ' && '.join([f"cp tmp1Warp.nii.gz {output.warp}", f"cp tmp1InverseWarp.nii.gz {output.invwarp}"]) @@ -907,8 +1000,8 @@ rule unfoldreg: label="{label}", **inputs.subj_wildcards ), -# shadow: -# "minimal" + shadow: + "minimal" threads: 16 resources: mem_mb=16000, @@ -979,8 +1072,10 @@ rule convert_unfoldreg_warp_from_itk_to_world: suffix="xfm.nii.gz", datatype="warps", desc="{desc}", - from_="{to}", #swapped around since surf needs opposite direction - to="{from}", + # from_="{to}", #swapped around since surf needs opposite direction + # to="{from}", + from_="{from}", + to="{to}", space="{space}", type_="itk", hemi="{hemi}", @@ -1007,6 +1102,30 @@ rule convert_unfoldreg_warp_from_itk_to_world: "wb_command -convert-warpfield -from-itk {input} -to-world {output}" +def get_unfold_ref(wildcards): + """ function to return either unfoldreg or unfold ref mesh, depending on whether + unfoldreg can be performed (based on atlas wildcards)""" + + if wildcards.label in config['atlas_files'][config['atlas']]['label_wildcards'] and config['no_unfolded_reg'] == False: + return bids( + root=work, + datatype="surf", + suffix="midthickness.surf.gii", + space="unfoldreg", + hemi="{hemi}", + label="{label}", + **inputs.subj_wildcards + ) + else: + return bids( + root=work, + datatype="surf", + suffix="midthickness.surf.gii", + space="unfold", + hemi="{hemi}", + label="{label}", + **inputs.subj_wildcards + ) rule warp_unfold_native_to_unfoldreg: @@ -1027,7 +1146,7 @@ rule warp_unfold_native_to_unfoldreg: datatype="warps", desc="SyN", from_="native", - to=config['atlas'][0], + to=config['atlas'], space="unfold", type_="surface", hemi="{hemi}", @@ -1057,49 +1176,33 @@ rule warp_unfold_native_to_unfoldreg: "wb_command -set-structure {output.surf_gii} {params.structure_type} -surface-type {params.surface_type}" " -surface-secondary-type {params.secondary_type}" - - -#now we have native vertices, warped to unfoldreg atlas space, so we can now use these to resample our meshes and metrics -rule resample_subfields_to_native_surf: +# --- resampling using the unfoldreg surface to standard densities (0p5mm, 1mm, 2mm, unfoldiso) +# +rule resample_atlas_subfields_to_std_density: + """ resamples subfields from a custom atlas mesh (e.g. from an atlas anatomical/native mesh + warped to unfold space) to a standard density""" input: - label_gii=bids( - root=work, - datatype="surf", - den=resample_from_density, - suffix="subfields.label.gii", - space="unfold", - hemi="{hemi}", - label="hipp", - atlas="{atlas}", - **inputs.subj_wildcards - ), - ref_unfold=os.path.join( #this uses the standard unfoldiso space - we should replace this with an atlas-specific gifti + atlas_dir=lambda wildcards: Path(download_dir) / "atlas" / wildcards.atlas, + ref_unfold=os.path.join( workflow.basedir, "..", "resources", - "unfold_template_hipp", - "tpl-avg_space-unfold_den-{density}_midthickness.surf.gii".format( - density=resample_from_density - ), - ), - native_unfold=bids( - root=work, - datatype="surf", - suffix="midthickness.surf.gii", - space="unfoldreg", - atlas="{atlas}", - hemi="{hemi}", - label="{label}", - **inputs.subj_wildcards + "unfold_template_{label}", + "tpl-avg_space-unfold_den-{density}_midthickness.surf.gii", ), + params: + label_gii = lambda wildcards, input: Path(input.atlas_dir) / config["atlas_files"][wildcards.atlas]['label_gii'].format(**wildcards), + atlas_unfold=lambda wildcards, input: Path(input.atlas_dir) / config["atlas_files"][wildcards.atlas]['surf_gii'].format(surf_type='midthickness',**wildcards), + output: label_gii=bids( - root=work, + root=root, datatype="surf", suffix="subfields.label.gii", - space="corobl", + space="{space}", + den="{density}", hemi="{hemi}", - label="hipp", + label="{label}", atlas="{atlas}", **inputs.subj_wildcards ), @@ -1109,10 +1212,12 @@ rule resample_subfields_to_native_surf: group: "subj" shell: - "wb_command -label-resample {input.label_gii} {input.ref_unfold} {input.native_unfold} BARYCENTRIC {output.label_gii} -bypass-sphere-check" + "wb_command -label-resample {params.label_gii} {params.atlas_unfold} {input.ref_unfold} BARYCENTRIC {output.label_gii} -bypass-sphere-check" + -rule resample_native_surf: + +rule resample_native_surf_to_std_density: input: native=bids( root=work, @@ -1130,22 +1235,14 @@ rule resample_native_surf: "unfold_template_{label}", "tpl-avg_space-unfold_den-{density}_midthickness.surf.gii", ), - native_unfold=bids( - root=work, - datatype="surf", - suffix="midthickness.surf.gii", - space="unfoldreg", - hemi="{hemi}", - label="{label}", - **inputs.subj_wildcards - ), + native_unfold=get_unfold_ref output: native_resampled=bids( root=work, datatype="surf", suffix="{surf_name,midthickness}.surf.gii", - space="{space}", + space="{space,unfoldreg|corobl}", den="{density}", hemi="{hemi}", label="{label}", @@ -1159,8 +1256,7 @@ rule resample_native_surf: shell: "wb_command -surface-resample {input.native} {input.native_unfold} {input.ref_unfold} BARYCENTRIC {output.native_resampled} -bypass-sphere-check" - -rule resample_native_metric: +rule resample_native_metric_to_std_density: input: native_metric=bids( root=work, @@ -1178,15 +1274,7 @@ rule resample_native_metric: "unfold_template_{label}", "tpl-avg_space-unfold_den-{density}_midthickness.surf.gii", ), - native_unfold=bids( - root=work, - datatype="surf", - suffix="midthickness.surf.gii", - space="unfoldreg", - hemi="{hemi}", - label="{label}", - **inputs.subj_wildcards - ), + native_unfold=get_unfold_ref output: metric_resampled=bids( @@ -1212,7 +1300,7 @@ rule cp_surf_to_root: native_resampled=bids( root=work, datatype="surf", - suffix="{surf_name,midthickness}.surf.gii", + suffix="{surf_name}.surf.gii", space="{space}", den="{density}", hemi="{hemi}", @@ -1232,31 +1320,80 @@ rule cp_surf_to_root: ), shell: "cp {input} {output}" -""" -rule cp_metric_to_root: - input: - metric_resampled=bids( + + +# --- resampling from atlasnative to native vertices +rule resample_atlas_subfields_to_native_surf: + input: + atlas_dir=lambda wildcards: Path(download_dir) / "atlas" / wildcards.atlas, + native_unfold=get_unfold_ref + params: + label_gii = lambda wildcards, input: Path(input.atlas_dir) / config["atlas_files"][wildcards.atlas]['label_gii'].format(**wildcards), + ref_unfold=lambda wildcards, input: Path(input.atlas_dir) / config["atlas_files"][wildcards.atlas]['surf_gii'].format(surf_type='midthickness',**wildcards), + + output: + label_gii=bids( root=work, datatype="surf", - suffix="{metricname}.{metrictype}.gii", - space="{space}", - den="{density}", - atlas=config['atlas'][0], + suffix="subfields.label.gii", + space="corobl", hemi="{hemi}", + label="hipp", + atlas="{atlas}", + **inputs.subj_wildcards + ), + container: + None + # config["singularity"]["autotop"] + group: + "subj" + shell: + "wb_command -label-resample {params.label_gii} {params.ref_unfold} {input.native_unfold} BARYCENTRIC {output.label_gii} -bypass-sphere-check" + +# --- rule for converting atlasnative subfields to unfold nii (e.g. analogous to unfoldiso standard space) +rule atlas_label_to_unfold_nii: + """converts metric .gii files to .nii - this is a band-aid fix since we make use of the volumetric + subfield labels in unfold space for painting the native volumetric dseg. Ie this uses the unfold volumetric (or effectively unfoldiso) + to perform mapping from atlas to subject. better approach would be to adjust the downstream volumetric dseg function to + make use of gifti labels instead.. + ALSO - this has issues (with ribbon mapping) since the most inner and outer coords are missing.. + -- maybe nearest-vertex works better? +""" + + input: + atlas_dir=lambda wildcards: Path(download_dir) / "atlas" / wildcards.atlas, + ref_nii=bids( + root=root, + space="unfold", label="{label}", - **inputs.subj_wildcards, + datatype="warps", + suffix="refvol.nii.gz", + **inputs.subj_wildcards ), + params: + label_gii = lambda wildcards, input: Path(input.atlas_dir) / config["atlas_files"][wildcards.atlas]['label_gii'].format(**wildcards), + inner_surf=lambda wildcards, input: Path(input.atlas_dir) / config["atlas_files"][wildcards.atlas]['surf_gii'].format(surf_type='inner',**wildcards), + outer_surf=lambda wildcards, input: Path(input.atlas_dir) / config["atlas_files"][wildcards.atlas]['surf_gii'].format(surf_type='outer',**wildcards), + midthickness_surf=lambda wildcards, input: Path(input.atlas_dir) / config["atlas_files"][wildcards.atlas]['surf_gii'].format(surf_type='midthickness',**wildcards), output: - metric_resampled=bids( - root=root, - datatype="surf", - suffix="{metricname}.{metrictype,shape|func}.gii", - space="{space}", - den="{density}", + label_nii=bids( + root=work, + datatype="anat", + suffix="subfields.nii.gz", + space="unfold", hemi="{hemi}", label="{label}", - **inputs.subj_wildcards, + atlas='{atlas}', + **inputs.subj_wildcards ), - shell: 'cp {input} {output}' + container: + config["singularity"]["autotop"] + group: + "subj" + shell: + "wb_command -label-to-volume-mapping {params.label_gii} {params.midthickness_surf} {input.ref_nii} {output.label_nii} " + " -nearest-vertex 1000" #use really large distance to ensure all voxels labelled +# " -ribbon-constrained {params.inner_surf} {params.outer_surf}" + + -""" diff --git a/hippunfold/workflow/rules/resample_final_to_crop_native.smk b/hippunfold/workflow/rules/resample_final_to_crop_native.smk index 7402c148..fcbd85ea 100644 --- a/hippunfold/workflow/rules/resample_final_to_crop_native.smk +++ b/hippunfold/workflow/rules/resample_final_to_crop_native.smk @@ -13,7 +13,7 @@ rule create_native_crop_ref: atlas="{atlas}", **inputs.subj_wildcards ), - atlas=config["atlas"][0], + atlas=config["atlas"], allow_missing=True, ), params: diff --git a/hippunfold/workflow/rules/subfields.smk b/hippunfold/workflow/rules/subfields.smk index 5af5a7e9..9b6845e5 100644 --- a/hippunfold/workflow/rules/subfields.smk +++ b/hippunfold/workflow/rules/subfields.smk @@ -1,54 +1,8 @@ - -rule resample_unfoldreg_subfields: +#TODO: replace this with a function that operates on the label gifti instead for native_warp workflow +# to avoid use of the unfold2native warps +rule label_subfields_from_vol_coords_corobl: + """ Label subfields using the volumetric coords and atlas subfield labels""" input: - label_nii=bids( - root=work, - datatype="anat", - suffix="subfields.nii.gz", - space="unfold", - hemi="{hemi}", - label="hipp", - atlas="{atlas}", - **inputs.subj_wildcards - ), - warp=bids( - root=work, - **inputs.subj_wildcards, - suffix="xfm.nii.gz", - datatype="warps", - desc="SyN", - from_="{atlas}", - to="subject", - space="unfold", - type_="itk", - hemi="{hemi}" - ), - output: - label_nii=bids( - root=root, - datatype="anat", - suffix="subfields.nii.gz", - space="unfold", - hemi="{hemi}", - label="hipp", - atlas="{atlas}", - **inputs.subj_wildcards - ), - container: - config["singularity"]["autotop"] - shadow: - "minimal" - group: - "subj" - shell: - "c3d {input.label_nii} -slice z 0:15 -oo tmp0%d.nii.gz && " - "for fn in $(ls tmp*.nii.gz); do antsApplyTransforms -d 2 -i $fn -r $fn -o $fn -n MultiLabel -t {input.warp}; done && " - "c3d tmp*.nii.gz -tile z -o recombined.nii.gz && " - "c3d {input.label_nii} recombined.nii.gz -copy-transform -o {output}" - - -def skip_unfoldreg_option_subfields(wildcards): - if config["no_unfolded_reg"]: label_nii = bids( root=work, datatype="anat", @@ -58,27 +12,7 @@ def skip_unfoldreg_option_subfields(wildcards): label="hipp", atlas="{atlas}", **inputs.subj_wildcards - ) - else: - label_nii = bids( - root=root, - datatype="anat", - suffix="subfields.nii.gz", - space="unfold", - hemi="{hemi}", - label="hipp", - atlas="{atlas}", - **inputs.subj_wildcards - ) - return label_nii - - -#TODO: replace this with a function that operates on the label gifti instead for native_warp workflow -# to avoid use of the unfold2native warps -rule label_subfields_from_vol_coords_corobl: - """ Label subfields using the volumetric coords and atlas subfield labels""" - input: - label_nii=skip_unfoldreg_option_subfields, + ), nii_ap=bids( root=work, datatype="coords", @@ -131,7 +65,7 @@ def get_tissue_atlas_remapping(wildcards): for label in mapping["tissue"].keys(): in_label = mapping["tissue"][label] - out_label = mapping[wildcards.atlas][label] + out_label = mapping.get(wildcards.atlas,mapping.get('default'))[label] remap.append(f"-threshold {in_label} {in_label} {out_label} 0 -popas {label}") diff --git a/hippunfold/workflow/rules/unfoldreg_v1.smk b/hippunfold/workflow/rules/unfoldreg_v1.smk index 4682c1dc..793d0b8d 100644 --- a/hippunfold/workflow/rules/unfoldreg_v1.smk +++ b/hippunfold/workflow/rules/unfoldreg_v1.smk @@ -858,5 +858,53 @@ rule cp_unfolded_noconstrain: shell: "cp {input} {output}" +#--- from subfields.smk + +rule resample_unfoldreg_subfields: + input: + label_nii=bids( + root=work, + datatype="anat", + suffix="subfields.nii.gz", + space="unfold", + hemi="{hemi}", + label="hipp", + atlas="{atlas}", + **inputs.subj_wildcards + ), + warp=bids( + root=work, + **inputs.subj_wildcards, + suffix="xfm.nii.gz", + datatype="warps", + desc="SyN", + from_="{atlas}", + to="subject", + space="unfold", + type_="itk", + hemi="{hemi}" + ), + output: + label_nii=bids( + root=root, + datatype="anat", + suffix="subfields.nii.gz", + space="unfold", + hemi="{hemi}", + label="hipp", + atlas="{atlas}", + **inputs.subj_wildcards + ), + container: + config["singularity"]["autotop"] + shadow: + "minimal" + group: + "subj" + shell: + "c3d {input.label_nii} -slice z 0:15 -oo tmp0%d.nii.gz && " + "for fn in $(ls tmp*.nii.gz); do antsApplyTransforms -d 2 -i $fn -r $fn -o $fn -n MultiLabel -t {input.warp}; done && " + "c3d tmp*.nii.gz -tile z -o recombined.nii.gz && " + "c3d {input.label_nii} recombined.nii.gz -copy-transform -o {output}" diff --git a/hippunfold/workflow/rules/warps.smk b/hippunfold/workflow/rules/warps.smk index e2e62854..6c02954d 100644 --- a/hippunfold/workflow/rules/warps.smk +++ b/hippunfold/workflow/rules/warps.smk @@ -553,233 +553,3 @@ rule create_warps_dentate: "../scripts/create_warps.py" -rule expand_unfolded_warps: - """unfolded space registration in 2D expanded to 3D""" - input: - warp2d=bids( - root=work, - **inputs.subj_wildcards, - suffix="xfm.nii.gz", - datatype="warps", - desc="SyN", - from_="{from}", - to="{to}", - space="unfold", - type_="itk", - hemi="{hemi}" - ), - unfold_phys_coords_nii=bids( - root=work, - space="unfold", - label="hipp", - datatype="coords", - suffix="refcoords.nii.gz", - **inputs.subj_wildcards - ), - output: - warp3d=bids( - root=work, - **inputs.subj_wildcards, - suffix="xfm.nii.gz", - datatype="warps", - desc="SyN3D", - from_="{from}", - to="{to}", - space="unfold", - type_="itk", - hemi="{hemi}" - ), - group: - "subj" - container: - config["singularity"]["autotop"] - script: - "../scripts/expand_2Dwarp.py" - - -def get_unfold2unfoldatlas(wildcards): - if config["no_unfolded_reg"]: - fn = [] - else: - fn = ( - bids( - root=work, - **inputs.subj_wildcards, - suffix="xfm.nii.gz", - datatype="warps", - desc="SyN3D", - from_="subject", - to=config["atlas"][0], - space="unfold", - type_="itk", - hemi="{hemi}" - ), - ) - return fn - - -rule compose_warps_native_to_unfold: - """ Compose warps from native to unfold """ - input: - unfold2unfoldatlas=get_unfold2unfoldatlas, - corobl2unfold=bids( - root=work, - datatype="warps", - **inputs.subj_wildcards, - label="{autotop}", - suffix="xfm.nii.gz", - hemi="{hemi}", - from_="corobl", - to="unfold", - mode="image" - ), - ref=bids( - root=root, - space="unfold", - label="{autotop}", - datatype="warps", - suffix="refvol.nii.gz", - **inputs.subj_wildcards - ), - native2corobl=bids( - root=work, - datatype="warps", - **inputs.subj_wildcards, - suffix="xfm.txt", - from_="{native_modality}", - to="corobl", - desc="affine", - type_="itk" - ), - output: - bids( - root=root, - datatype="warps", - **inputs.subj_wildcards, - label="{autotop}", - suffix="xfm.nii.gz", - hemi="{hemi}", - from_="{native_modality}", - to="unfold", - mode="image" - ), - log: - bids( - root="logs", - **inputs.subj_wildcards, - label="{autotop}", - suffix="composexfm.txt", - hemi="{hemi}", - from_="{native_modality}", - to="unfold", - mode="image" - ), - container: - config["singularity"]["autotop"] - group: - "subj" - shell: - "ComposeMultiTransform 3 {output} -R {input.ref} {input.unfold2unfoldatlas} {input.corobl2unfold} {input.native2corobl} &> {log}" - - -def get_unfoldatlas2unfold(wildcards): - if config["no_unfolded_reg"]: - fn = [] - else: - fn = ( - bids( - root=work, - **inputs.subj_wildcards, - suffix="xfm.nii.gz", - datatype="warps", - desc="SyN3D", - from_=config["atlas"][0], - to="subject", - space="unfold", - type_="itk", - hemi="{hemi}", - ), - ) - return fn - - -def get_cmd_compose_warps_unfold_to_crop_native(wildcards, input, output): - if config["no_unfolded_reg"]: - cmd = f"antsApplyTransforms -o [{output.unfold2cropnative},1] -r {input.ref} -t [{input.native2corobl},1] -t {input.unfold2corobl} -i {input.unfold_ref} -v" - else: - cmd = f"antsApplyTransforms -o [{output.unfold2cropnative},1] -r {input.ref} -t [{input.native2corobl},1] -t {input.unfold2corobl} -t {input.unfoldatlas2unfold} -i {input.unfold_ref} -v" - return cmd - - -rule compose_warps_unfold_to_crop_native: - """ Compose warps from unfold to crop native """ - input: - unfoldatlas2unfold=get_unfoldatlas2unfold, - unfold2corobl=bids( - root=work, - datatype="warps", - **inputs.subj_wildcards, - label="{autotop}", - suffix="xfm.nii.gz", - hemi="{hemi}", - from_="unfold", - to="corobl", - mode="image" - ), - ref=bids( - root=work, - datatype="warps", - suffix="cropref.nii.gz", - space="{native_modality}", - hemi="{hemi}", - **inputs.subj_wildcards - ), - native2corobl=bids( - root=work, - datatype="warps", - **inputs.subj_wildcards, - suffix="xfm.txt", - from_="{native_modality}", - to="corobl", - desc="affine", - type_="itk" - ), - unfold_ref=bids( - root=root, - space="unfold", - label="{autotop}", - datatype="warps", - suffix="refvol.nii.gz", - **inputs.subj_wildcards - ), - output: - unfold2cropnative=bids( - root=root, - datatype="warps", - **inputs.subj_wildcards, - label="{autotop}", - suffix="xfm.nii.gz", - hemi="{hemi}", - from_="unfold", - to="{native_modality}", - mode="image" - ), - log: - bids( - root="logs", - **inputs.subj_wildcards, - label="{autotop}", - suffix="composexfm.txt", - hemi="{hemi}", - from_="unfold", - to="{native_modality}", - mode="image" - ), - params: - cmd=get_cmd_compose_warps_unfold_to_crop_native, - container: - config["singularity"]["autotop"] - group: - "subj" - shell: - "{params.cmd} &> {log}"