diff --git a/hippunfold/config/snakebids.yml b/hippunfold/config/snakebids.yml index 7d842739..e5a654a8 100644 --- a/hippunfold/config/snakebids.yml +++ b/hippunfold/config/snakebids.yml @@ -122,12 +122,6 @@ parse_args: --inject_template: choices: - 'upenn' - - 'dHCP' - - 'MBMv2' - - 'MBMv3' - - 'CIVM' - - 'ABAv3' - - 'bigbrain' default: 'upenn' help: 'Set the template to use for shape injection. (default: %(default)s)' @@ -189,6 +183,16 @@ parse_args: help: 'Resample the manual segmentation (for --modality dsegtissue), keeping the bounding box the same, but changing the number of voxels in the image. This can be specified as voxels ("300x300x300"), or a percentage ("20%%") or percentage in each direction ("20x20x40%%"). Note: the segmentation will always be subsequently cropped around the segmentation, with a padding of 5 voxels. (default: %(default)s)' default: null + --laminar_coords_res_hipp: + help: 'Set the resolution that the laminar (IO) coords will be computed with for hippocampus. This is implemented by resampling the segmentation to this resolution prior to generating the coords. (default: %(default)s)' + default: '0.3x0.3x0.3mm' + dest: "laminar_coords_res.hipp" + + --laminar_coords_res_dentate: + help: 'Set the resolution that the laminar (IO) coords will be computed with for the dentate gyrus. This is implemented by resampling the segmentation to this resolution prior to generating the coords. (default: %(default)s)' + default: '0.1x0.1x0.1mm' + dest: "laminar_coords_res.dentate" + --atlas: choices: @@ -426,10 +430,18 @@ template_files: upenn: T1w: tpl-upenn_desc-hipptissue_dseg.nii.gz dseg: tpl-upenn_desc-hipptissue_dseg.nii.gz + dseg_dentate: tpl-upenn_desc-hippdentate_dseg.nii.gz # hack: separate files for DG coords: tpl-upenn_dir-{dir}_label-{label}_coords.nii.gz hemi_wildcards: - R + bbhist: + T1w: sub-bbhist_hemi-{hemi}_100um.nii.gz + dseg: sub-bbhist_hemi-{hemi}_100um_dseg.nii.gz # aka dsegtissue + hemi_wildcards: + - L + - R + ABAv3: T1w: tpl-ABAv3_T1w.nii.gz xfm_corobl: tpl-ABAv3_from-native_to-corobl_type-itk_affine.txt @@ -475,15 +487,15 @@ atlas_files: lut: desc-subfields_dseg.tsv surf_gii: sub-bbhist_hemi-{hemi}_space-unfold_label-{label}_{surf_type}.surf.gii metric_gii: sub-bbhist_hemi-{hemi}_space-corobl_label-{label}_{metric}.shape.gii - label_gii: sub-bbhist_hemi-{hemi}_space-corobl_label-{label}_desc-manual_subfields.label.gii + label_gii: sub-bbhist_hemi-{hemi}_space-corobl_label-{label}_subfields.label.gii label_wildcards: - hipp + - dentate hemi_wildcards: - L - R metric_wildcards: - thickness - - gyrification - curvature bigbrain: @@ -570,6 +582,7 @@ resource_urls: CIVM: 'files.ca-1.osf.io/v1/resources/v8acf/providers/osfstorage/65395bf62827451220b86e24/?zip=' ABAv3: 'files.ca-1.osf.io/v1/resources/v8acf/providers/osfstorage/6668855b6b6c8e2cc704ca97/?zip=' bigbrain: 'files.ca-1.osf.io/v1/resources/v8acf/providers/osfstorage/666b1bc765e1de5972893e59/?zip=' + bbhist: 'placeholder.url' #to get hash, see https://github.com/CenterForOpenScience/osf.io/issues/8256#issuecomment-379833911 @@ -710,6 +723,14 @@ laplace_labels: - 9 - 10 - 0 + + +tight_crop_labels: #defines what labels to use for tight cropping (for upsampled ref) + hipp: + - 1 + dentate: + - 8 + output_spaces: - native @@ -737,3 +758,8 @@ workdir: null use_template_seg: False template_seg_smoothing_factor: 2 resample_dsegtissue: null +laminar_coords_res: + hipp: 0.3x0.3x0.3mm + dentate: 0.1x0.1x0.1mm +resample_dseg_for_templatereg: 0.3x0.3x0.3mm + diff --git a/hippunfold/workflow/rules/common.smk b/hippunfold/workflow/rules/common.smk index 0272df6a..cbfcca80 100644 --- a/hippunfold/workflow/rules/common.smk +++ b/hippunfold/workflow/rules/common.smk @@ -83,6 +83,7 @@ def get_final_subfields(): space="{space}", hemi="{hemi}", atlas="{atlas}", + label="hipp", **inputs.subj_wildcards, ), hemi=config["hemi"], @@ -331,3 +332,26 @@ def get_create_template_output(): ) ) return files + + +def get_input_for_shape_inject(wildcards): + if config["modality"] == "dsegtissue": + seg = bids( + root=work, + datatype="anat", + **inputs.subj_wildcards, + suffix="dseg.nii.gz", + space="corobl", + hemi="{hemi}", + ).format(**wildcards) + else: + seg = bids( + root=work, + datatype="anat", + **inputs.subj_wildcards, + suffix="dseg.nii.gz", + desc="nnunet", + space="corobl", + hemi="{hemi}", + ).format(**wildcards) + return seg diff --git a/hippunfold/workflow/rules/coords.smk b/hippunfold/workflow/rules/coords.smk index e93711d9..e2b25e23 100644 --- a/hippunfold/workflow/rules/coords.smk +++ b/hippunfold/workflow/rules/coords.smk @@ -13,6 +13,7 @@ def get_labels_for_laplace(wildcards): desc="postproc", space="corobl", hemi="{hemi}", + label="{label}", ).format(**wildcards) return seg @@ -185,6 +186,36 @@ rule get_nan_mask: "c3d {input} -background -1 -retain-labels {params} -binarize {output}" +rule create_upsampled_coords_ref: + input: + seg=get_input_for_shape_inject, + params: + tight_crop_labels=lambda wildcards: config["tight_crop_labels"][wildcards.label], + resample_res=lambda wildcards: config["laminar_coords_res"][wildcards.label], + trim_padding="3mm", + output: + upsampled_ref=bids( + root=work, + datatype="anat", + **inputs.subj_wildcards, + suffix="ref.nii.gz", + desc="resampled", + space="corobl", + label="{label}", + hemi="{hemi}", + ), + container: + config["singularity"]["autotop"] + conda: + "../envs/c3d.yaml" + group: + "subj" + container: + config["singularity"]["autotop"] + shell: + "c3d {input} -retain-labels {params.tight_crop_labels} -trim {params.trim_padding} -resample-mm {params.resample_res} -o {output}" + + rule prep_dseg_for_laynii: input: dseg_tissue=get_labels_for_laplace, diff --git a/hippunfold/workflow/rules/download.smk b/hippunfold/workflow/rules/download.smk index f1f83e7d..1691bdac 100644 --- a/hippunfold/workflow/rules/download.smk +++ b/hippunfold/workflow/rules/download.smk @@ -87,6 +87,43 @@ rule import_template_dseg: "{params.copy_or_flip_cmd} {output.template_seg}" +rule import_template_dseg_dentate: + input: + template_dir=Path(download_dir) / "template" / config["inject_template"], + params: + template_seg=lambda wildcards: Path(download_dir) + / "template" + / config["inject_template"] + / config["template_files"][config["inject_template"]]["dseg_dentate"].format( + **wildcards + ), + copy_or_flip_cmd=lambda wildcards: copy_or_flip( + wildcards, + Path(download_dir) + / "template" + / config["inject_template"] + / config["template_files"][config["inject_template"]][ + "dseg_dentate" + ].format(**wildcards), + ), + output: + template_seg=bids( + root=work, + datatype="anat", + space="template", + **inputs.subj_wildcards, + desc="dentatetissue", + hemi="{hemi}", + suffix="dseg.nii.gz", + ), + group: + "subj" + container: + config["singularity"]["autotop"] + shell: + "{params.copy_or_flip_cmd} {output.template_seg}" + + rule import_template_coords: input: template_dir=Path(download_dir) / "template" / config["inject_template"], diff --git a/hippunfold/workflow/rules/native_surf.smk b/hippunfold/workflow/rules/native_surf.smk index 3e691a9d..6d9471ae 100644 --- a/hippunfold/workflow/rules/native_surf.smk +++ b/hippunfold/workflow/rules/native_surf.smk @@ -5,6 +5,7 @@ surf_thresholds = {"inner": 0, "outer": 1, "midthickness": 0.5} + unfoldreg_method = "greedy" # choices: ["greedy","SyN"] unfoldreg_padding = "64x64x0vox" @@ -66,7 +67,7 @@ rule gen_native_mesh: params: threshold=lambda wildcards: surf_thresholds[wildcards.surfname], decimate_opts={ - "reduction": 0.9, + "reduction": 0.7, "feature_angle": 25, "preserve_topology": True, }, @@ -292,9 +293,12 @@ rule laplace_beltrami: **inputs.subj_wildcards, ), params: - min_dist_threshold=0.3, - max_dist_threshold=1, - min_terminal_vertices=5, + min_dist_percentile=1, + max_dist_percentile=10, + min_terminal_vertices=lambda wildcards: 5 if wildcards.dir == "AP" else 100, #TODO, instead of # of vertices, we should compute the total length of the segment + threshold_method=lambda wildcards: ( + "percentile" if wildcards.dir == "AP" else "firstminima" + ), output: coords=bids( root=work, @@ -692,12 +696,15 @@ rule warp_midthickness_to_inout: config["singularity"]["autotop"] conda: "../envs/workbench.yaml" + shadow: + "minimal" group: "subj" shell: - "wb_command -surface-apply-warpfield {input.surf_gii} {input.warp} {output.surf_gii} && " - "wb_command -set-structure {output.surf_gii} {params.structure_type} -surface-type {params.surface_type}" - " -surface-secondary-type {params.secondary_type}" + "wb_command -volume-to-surface-mapping {input.warp} {input.surf_gii} warp.shape.gii -trilinear && " + "wb_command -surface-coordinates-to-metric {input.surf_gii} coords.shape.gii && " + "wb_command -metric-math 'COORDS + WARP' warpedcoords.shape.gii -var COORDS coords.shape.gii -var WARP warp.shape.gii && " + "wb_command -surface-set-coordinates {input.surf_gii} warpedcoords.shape.gii {output.surf_gii}" # --- affine transforming anatomical surfaces from corobl to other (T1w, T2w) spaces diff --git a/hippunfold/workflow/rules/qc.smk b/hippunfold/workflow/rules/qc.smk index b284f2d5..0364f296 100644 --- a/hippunfold/workflow/rules/qc.smk +++ b/hippunfold/workflow/rules/qc.smk @@ -51,6 +51,7 @@ rule get_subfield_vols_subj: space="{crop_ref_spaces}", desc="subfields", atlas="{atlas}", + label="hipp", suffix="dseg.nii.gz", ), hemi=config["hemi"], @@ -173,6 +174,7 @@ rule qc_subfield: space="{space}", hemi="{hemi}", atlas="{atlas}", + label="hipp", **inputs.subj_wildcards, ), output: diff --git a/hippunfold/workflow/rules/resample_final_to_crop_native.smk b/hippunfold/workflow/rules/resample_final_to_crop_native.smk index b650a89f..5055dddf 100644 --- a/hippunfold/workflow/rules/resample_final_to_crop_native.smk +++ b/hippunfold/workflow/rules/resample_final_to_crop_native.smk @@ -11,6 +11,7 @@ rule create_native_crop_ref: space="{native_modality}", hemi="{hemi}", atlas="{atlas}", + label="hipp", **inputs.subj_wildcards, ), atlas=config["atlas"], @@ -98,6 +99,7 @@ rule resample_postproc_native_crop: desc="postproc", space="corobl", hemi="{hemi}", + label=config["autotop_labels"][-1], ), xfm=bids( root=work, @@ -148,6 +150,7 @@ rule resample_subfields_native_crop: space="corobl", hemi="{hemi}", atlas="{atlas}", + label="{label}", **inputs.subj_wildcards, ), xfm=bids( @@ -177,6 +180,7 @@ rule resample_subfields_native_crop: space="crop{native_modality}", hemi="{hemi}", atlas="{atlas}", + label="{label,hipp}", **inputs.subj_wildcards, ), container: diff --git a/hippunfold/workflow/rules/shape_inject.smk b/hippunfold/workflow/rules/shape_inject.smk index 349302ad..5ca067e7 100644 --- a/hippunfold/workflow/rules/shape_inject.smk +++ b/hippunfold/workflow/rules/shape_inject.smk @@ -2,29 +2,6 @@ # by flipping it -def get_input_for_shape_inject(wildcards): - if config["modality"] == "dsegtissue": - seg = bids( - root=work, - datatype="anat", - **inputs.subj_wildcards, - suffix="dseg.nii.gz", - space="corobl", - hemi="{hemi}", - ).format(**wildcards) - else: - seg = bids( - root=work, - datatype="anat", - **inputs.subj_wildcards, - suffix="dseg.nii.gz", - desc="nnunet", - space="corobl", - hemi="{hemi}", - ).format(**wildcards) - return seg - - def get_input_splitseg_for_shape_inject(wildcards): if config["modality"] == "dsegtissue": seg = bids( @@ -112,7 +89,7 @@ def get_inject_scaling_opt(wildcards): return f"-s {gradient_sigma}vox {warp_sigma}vox" -rule template_shape_reg: +rule resample_template_dseg_tissue_for_reg: input: template_seg=bids( root=work, @@ -121,12 +98,48 @@ rule template_shape_reg: **inputs.subj_wildcards, desc="hipptissue", hemi="{hemi}", + suffix="dseg.nii.gz", + ), + params: + resample_cmd="-resample-mm {res}".format( + res=config["resample_dseg_for_templatereg"] + ), + crop_cmd="-trim 5vox", #leave 5 voxel padding + output: + template_seg=bids( + root=work, + datatype="anat", + space="template", + **inputs.subj_wildcards, + desc="hipptissueresampled", + hemi="{hemi}", + suffix="dseg.nii.gz", + ), + container: + config["singularity"]["autotop"] + conda: + "../envs/c3d.yaml" + group: + "subj" + shell: + "c3d {input} -int 0 {params.resample_cmd} {params.crop_cmd} -o {output}" + + +rule template_shape_reg: + input: + template_seg=bids( + root=work, + datatype="anat", + space="template", + **inputs.subj_wildcards, + desc="hipptissueresampled", + hemi="{hemi}", suffix="dsegsplit", ), subject_seg=get_input_splitseg_for_shape_inject, params: general_opts="-d 3 -m SSD", - affine_opts="-moments 2", + affine_opts="-moments 2 -det 1", greedy_opts=get_inject_scaling_opt, img_pairs=get_image_pairs, output: @@ -173,7 +186,9 @@ rule template_shape_reg: "greedy -threads {threads} {params.general_opts} {params.greedy_opts} {params.img_pairs} -it {output.matrix} -o {output.warp} &>> {log}" -rule template_shape_inject: +rule dilate_dentate_pd_src_sink: + """ The PD src/sink labels can disappear after label propagation + as they are very small. This dilates them into relative background labels""" input: template_seg=bids( root=work, @@ -184,7 +199,53 @@ rule template_shape_inject: hemi="{hemi}", suffix="dseg.nii.gz", ), - subject_seg=get_input_for_shape_inject, + params: + src_label=config["laplace_labels"]["dentate"]["PD"]["src"][0], + sink_label=config["laplace_labels"]["dentate"]["PD"]["sink"][0], + src_bg=2, + sink_bg=10, + struc_elem_size=3, + output: + template_seg=bids( + root=work, + datatype="anat", + space="template", + **inputs.subj_wildcards, + desc="hipptissuedilated", + hemi="{hemi}", + suffix="dseg.nii.gz", + ), + group: + "subj" + container: + config["singularity"]["autotop"] + conda: + "../envs/neurovis.yaml" + script: + "../scripts/dilate_dentate_pd_src_sink.py" + + +rule template_shape_inject: + input: + template_seg=bids( + root=work, + datatype="anat", + space="template", + **inputs.subj_wildcards, + desc="{label}tissue", + hemi="{hemi}", + suffix="dseg.nii.gz", + ), + upsampled_ref=bids( + root=work, + datatype="anat", + **inputs.subj_wildcards, + suffix="ref.nii.gz", + desc="resampled", + label="{label}", + space="corobl", + hemi="{hemi}", + ), matrix=bids( root=work, **inputs.subj_wildcards, @@ -209,7 +270,7 @@ rule template_shape_inject: hemi="{hemi}", ), params: - interp_opt="-ri LABEL 0.2vox", + interp_opt="-ri LABEL 0.1mm", # smoothing sigma = 100micron output: inject_seg=bids( root=work, @@ -219,6 +280,7 @@ rule template_shape_inject: desc="inject", space="corobl", hemi="{hemi}", + label="{label}", ), log: bids( @@ -226,6 +288,7 @@ rule template_shape_inject: **inputs.subj_wildcards, suffix="templateshapeinject.txt", hemi="{hemi}", + label="{label}", ), group: "subj" @@ -235,12 +298,22 @@ rule template_shape_inject: "../envs/greedy.yaml" threads: 8 shell: - "greedy -d 3 -threads {threads} {params.interp_opt} -rf {input.subject_seg} -rm {input.template_seg} {output.inject_seg} -r {input.warp} {input.matrix} &> {log}" + "greedy -d 3 -threads {threads} {params.interp_opt} -rf {input.upsampled_ref} -rm {input.template_seg} {output.inject_seg} -r {input.warp} {input.matrix} &> {log}" rule inject_init_laplace_coords: + """ TODO: this may not be needed anymore """ input: - subject_seg=get_input_for_shape_inject, + upsampled_ref=bids( + root=work, + datatype="anat", + **inputs.subj_wildcards, + suffix="ref.nii.gz", + desc="resampled", + space="corobl", + hemi="{hemi}", + label="{label}", + ), matrix=bids( root=work, **inputs.subj_wildcards, @@ -276,7 +349,7 @@ rule inject_init_laplace_coords: hemi="{hemi}", ), params: - interp_opt="-ri NN", + interp_opt="-ri LIN", output: init_coords=bids( root=work, @@ -307,10 +380,15 @@ rule inject_init_laplace_coords: "../envs/greedy.yaml" threads: 8 shell: - "greedy -d 3 -threads {threads} {params.interp_opt} -rf {input.subject_seg} -rm {input.coords} {output.init_coords} -r {input.warp} {input.matrix} &> {log}" + "greedy -d 3 -threads {threads} {params.interp_opt} -rf {input.upsampled_ref} -rm {input.coords} {output.init_coords} -r {input.warp} {input.matrix} &> {log}" rule reinsert_subject_labels: + """ c3d command to: + 1) get the labels to retain + 2) reslice to injected seg + 3) set injected seg to zero where retained labels are + 4) and add this to retained labels""" input: inject_seg=bids( root=work, @@ -320,6 +398,7 @@ rule reinsert_subject_labels: desc="inject", space="corobl", hemi="{hemi}", + label="{label}", ), subject_seg=get_input_for_shape_inject, params: @@ -335,6 +414,7 @@ rule reinsert_subject_labels: desc="postproc", space="corobl", hemi="{hemi}", + label="{label}", ), group: "subj" @@ -343,4 +423,7 @@ rule reinsert_subject_labels: conda: "../envs/c3d.yaml" shell: - "c3d {input.subject_seg} -retain-labels {params.labels} -popas LBL -push LBL -threshold 0 0 1 0 {input.inject_seg} -multiply -push LBL -add -o {output.postproc_seg}" + "c3d {input.subject_seg} -retain-labels {params.labels} -popas LBL " + " -int 0 {input.inject_seg} -as SEG -push LBL -reslice-identity -popas LBL_RESLICE " + "-push LBL_RESLICE -threshold 0 0 1 0 -push SEG -multiply " + "-push LBL_RESLICE -add -o {output.postproc_seg}" diff --git a/hippunfold/workflow/rules/subfields.smk b/hippunfold/workflow/rules/subfields.smk index 54b092a0..97890e5c 100644 --- a/hippunfold/workflow/rules/subfields.smk +++ b/hippunfold/workflow/rules/subfields.smk @@ -50,7 +50,7 @@ rule subfields_to_label_gifti: suffix="subfields.label.gii", space="corobl", hemi="{hemi}", - label="hipp", + label="{label,hipp}", **inputs.subj_wildcards, ), conda: @@ -157,7 +157,7 @@ rule combine_tissue_subfield_labels_corobl: space="corobl", hemi="{hemi}", atlas="{atlas}", - label="hipp", + label="{label}", **inputs.subj_wildcards, ), params: @@ -170,6 +170,7 @@ rule combine_tissue_subfield_labels_corobl: suffix="dseg.nii.gz", space="corobl", hemi="{hemi}", + label="{label,hipp}", atlas="{atlas}", **inputs.subj_wildcards, ), @@ -194,6 +195,7 @@ rule resample_subfields_to_native: space="corobl", hemi="{hemi}", atlas="{atlas}", + label="{label}", **inputs.subj_wildcards, ), xfm=bids( @@ -222,6 +224,7 @@ rule resample_subfields_to_native: space="{native_modality,T1w|T2w}", hemi="{hemi}", atlas="{atlas}", + label="{label,hipp}", **inputs.subj_wildcards, ), container: @@ -246,6 +249,7 @@ rule resample_postproc_to_native: desc="postproc", space="corobl", hemi="{hemi}", + label=config["autotop_labels"][-1], ), xfm=bids( root=work, diff --git a/hippunfold/workflow/rules/templateseg.smk b/hippunfold/workflow/rules/templateseg.smk index a02da2ee..9a9ccfd5 100644 --- a/hippunfold/workflow/rules/templateseg.smk +++ b/hippunfold/workflow/rules/templateseg.smk @@ -66,14 +66,15 @@ rule template_reg: rule warp_template_dseg: input: - ref=bids( + upsampled_ref=bids( root=work, datatype="anat", **inputs.subj_wildcards, - suffix=f"{config['modality']}.nii.gz", + suffix="ref.nii.gz", + desc="resampled", space="corobl", - desc="preproc", hemi="{hemi}", + label="{label}", ), warp=bids( root=work, @@ -106,6 +107,7 @@ rule warp_template_dseg: desc="postproc", space="corobl", hemi="{hemi}", + label="{label}", ), group: "subj" diff --git a/hippunfold/workflow/scripts/dilate_dentate_pd_src_sink.py b/hippunfold/workflow/scripts/dilate_dentate_pd_src_sink.py new file mode 100644 index 00000000..e1b9203f --- /dev/null +++ b/hippunfold/workflow/scripts/dilate_dentate_pd_src_sink.py @@ -0,0 +1,45 @@ +import numpy as np +import nibabel as nib +from scipy.ndimage import binary_dilation, generate_binary_structure + + +def selective_dilation( + input_nifti, output_nifti, src, sink, src_bg, sink_bg, structure_size=3 +): + # Load image + img = nib.load(input_nifti) + data = img.get_fdata().astype(np.int32) + + # Create structuring element + structure = generate_binary_structure(3, 1) + + # Dilation: src -> src_bg + src_mask = data == src + src_bg_mask = data == src_bg + src_dilated = binary_dilation(src_mask, structure=structure) & src_bg_mask + data[src_dilated] = src + + # Dilation: sink -> sink_bg + sink_mask = data == sink + sink_bg_mask = data == sink_bg + sink_dilated = binary_dilation(sink_mask, structure=structure) & sink_bg_mask + data[sink_dilated] = sink + + # Save the modified image + new_img = nib.Nifti1Image(data, img.affine, img.header) + nib.save(new_img, output_nifti) + + print(f"Output saved to {output_nifti}") + + +input_nifti = snakemake.input.template_seg +output_nifti = snakemake.output.template_seg +src = snakemake.params.src_label +sink = snakemake.params.sink_label +src_bg = snakemake.params.src_bg +sink_bg = snakemake.params.sink_bg +structure_size = snakemake.params.struc_elem_size + +selective_dilation( + input_nifti, output_nifti, src, sink, src_bg, sink_bg, structure_size +) diff --git a/hippunfold/workflow/scripts/laplace_beltrami.py b/hippunfold/workflow/scripts/laplace_beltrami.py index e9cfabc8..ed953c3b 100644 --- a/hippunfold/workflow/scripts/laplace_beltrami.py +++ b/hippunfold/workflow/scripts/laplace_beltrami.py @@ -6,8 +6,107 @@ log_file = snakemake.log[0] if snakemake.log else None logger = setup_logger(log_file) +import numpy as np + +import numpy as np +from scipy.signal import argrelextrema + + +def get_terminal_indices_firstminima( + sdt, min_vertices, boundary_mask, bins=100, smoothing_window=5 +): + """ + Gets the terminal (src/sink) vertex indices by determining an adaptive threshold + using the first local minimum of the histogram of `sdt` values. + + Parameters: + - sdt: Signed distance transform array. + - min_vertices: The minimum number of vertices required. + - boundary_mask: Boolean or binary mask indicating boundary regions. + - bins: Number of bins to use in the histogram (default: 100). + - smoothing_window: Window size for moving average smoothing (default: 5). + + Returns: + - indices: List of terminal vertex indices. -def get_terminal_indices(sdt, min_dist, max_dist, min_vertices, boundary_mask): + Raises: + - ValueError: If the minimum number of vertices cannot be found. + """ + + # Extract SDT values within the boundary mask + sdt_values = sdt[boundary_mask == 1] + + # Compute histogram + hist, bin_edges = np.histogram(sdt_values, bins=bins, density=True) + bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2 + + # Smooth the histogram using a simple moving average + smoothed_hist = np.convolve( + hist, np.ones(smoothing_window) / smoothing_window, mode="same" + ) + + # Find local minima + minima_indices = argrelextrema(smoothed_hist, np.less)[0] + + if len(minima_indices) == 0: + raise ValueError("No local minima found in the histogram.") + + # Select the first local minimum after the first peak + first_minimum_bin = bin_centers[minima_indices[0]] + + # Select indices where SDT is below this threshold + indices = np.where((sdt < first_minimum_bin) & (boundary_mask == 1))[0].tolist() + + if len(indices) >= min_vertices: + return indices + + raise ValueError( + f"Unable to find minimum of {min_vertices} vertices on boundary within the first local minimum threshold." + ) + + +def get_terminal_indices_percentile( + sdt, min_percentile, max_percentile, min_vertices, boundary_mask +): + """ + Gets the terminal (src/sink) vertex indices by sweeping a percentile-based threshold + of the signed distance transform (sdt), ensuring at least `min_vertices` are selected. + + Instead of a fixed distance range, this function dynamically determines the threshold + by scanning from `min_percentile` to `max_percentile`. + + Parameters: + - sdt: Signed distance transform array. + - min_percentile: Starting percentile for thresholding (0-100). + - max_percentile: Maximum percentile for thresholding (0-100). + - min_vertices: The minimum number of vertices required. + - boundary_mask: Boolean or binary mask indicating boundary regions. + + Returns: + - indices: List of terminal vertex indices. + + Raises: + - ValueError: If the minimum number of vertices cannot be found. + """ + + for percentile in np.arange(min_percentile, max_percentile, 0.5): + dist_threshold = np.percentile(sdt[boundary_mask == 1], percentile) + indices = np.where((sdt < dist_threshold) & (boundary_mask == 1))[0].tolist() + + if len(indices) >= min_vertices: + logger.info( + f"Using {percentile}-th percentile to obtain sdt threshold of {dist_threshold}, with {len(indices)} vertices" + ) + return indices + + raise ValueError( + f"Unable to find minimum of {min_vertices} vertices on boundary within the {max_percentile}th percentile of distances" + ) + + +def get_terminal_indices_threshold( + sdt, min_dist, max_dist, min_vertices, boundary_mask +): """ Gets the terminal (src/sink) vertex indices based on distance to the src/sink mask, a boundary mask, and a minumum number of vertices. The distance from the mask is @@ -116,25 +215,40 @@ def solve_laplace_beltrami_open_mesh(vertices, faces, boundary_conditions=None): src_sdt = nib.load(snakemake.input.src_sdt).agg_data() sink_sdt = nib.load(snakemake.input.sink_sdt).agg_data() -src_indices = get_terminal_indices( - src_sdt, - snakemake.params.min_dist_threshold, - snakemake.params.max_dist_threshold, - snakemake.params.min_terminal_vertices, - boundary_mask, -) +if snakemake.params.threshold_method == "percentile": + src_indices = get_terminal_indices_percentile( + src_sdt, + snakemake.params.min_dist_percentile, + snakemake.params.max_dist_percentile, + snakemake.params.min_terminal_vertices, + boundary_mask, + ) + sink_indices = get_terminal_indices_percentile( + sink_sdt, + snakemake.params.min_dist_percentile, + snakemake.params.max_dist_percentile, + snakemake.params.min_terminal_vertices, + boundary_mask, + ) -logger.info(f"# of src boundary vertices: {len(src_indices)}") -sink_indices = get_terminal_indices( - sink_sdt, - snakemake.params.min_dist_threshold, - snakemake.params.max_dist_threshold, - snakemake.params.min_terminal_vertices, - boundary_mask, -) +elif snakemake.params.threshold_method == "firstminima": + src_indices = get_terminal_indices_firstminima( + src_sdt, + snakemake.params.min_terminal_vertices, + boundary_mask, + ) + sink_indices = get_terminal_indices_firstminima( + sink_sdt, + snakemake.params.min_terminal_vertices, + boundary_mask, + ) + + +logger.info(f"# of src boundary vertices: {len(src_indices)}") logger.info(f"# of sink boundary vertices: {len(sink_indices)}") + src_vals = [0 for i in range(len(src_indices))] sink_vals = [1 for i in range(len(sink_indices))]