From f4b918304b295ed243d9476f4842d6468d2d08c9 Mon Sep 17 00:00:00 2001 From: Ali Khan Date: Wed, 1 Jan 2025 23:59:32 -0500 Subject: [PATCH] fix for isosurface generation, and cleanup major fix for isosurfaces as was using wrong labels for nan (all background) when I should instead be useing the src/sink for the other coords (AP/PD). Also sets the sink label to 1.1 to eliminate the false boundary, and then also use nans to ignore the false boundary thats created where src and sink meet.. THis now produces hipp surfaces with no holes! --- hippunfold/workflow/rules/native_surf.smk | 155 ++++++++++-------- hippunfold/workflow/scripts/gen_isosurface.py | 64 ++++++-- .../prep_dentate_coords_for_meshing.py | 39 ----- .../scripts/prep_hipp_coords_for_meshing.py | 57 ------- 4 files changed, 141 insertions(+), 174 deletions(-) delete mode 100644 hippunfold/workflow/scripts/prep_dentate_coords_for_meshing.py delete mode 100644 hippunfold/workflow/scripts/prep_hipp_coords_for_meshing.py diff --git a/hippunfold/workflow/rules/native_surf.smk b/hippunfold/workflow/rules/native_surf.smk index 1b781cb6..1aa92954 100644 --- a/hippunfold/workflow/rules/native_surf.smk +++ b/hippunfold/workflow/rules/native_surf.smk @@ -12,6 +12,23 @@ gm_labels = { "dentate": config["laplace_labels"]["PD"]["sink"], } +# appends the coords with these regions set to +1.1 for the meshing +sink_labels = {"hipp": config["laplace_labels"]["IO"]["sink"], "dentate": [2]} + +# appends the coords with these regions set to +1.1 for the meshing +src_labels = {"hipp": config["laplace_labels"]["IO"]["src"], "dentate": [1]} + +# sets these to nan in the coords for the meshing +nan_labels = { + "hipp": config["laplace_labels"]["AP"]["sink"] + + config["laplace_labels"]["AP"]["src"] + + config["laplace_labels"]["PD"]["sink"] + + config["laplace_labels"]["PD"]["src"], + "dentate": [ + 0 + ], # TODO: this requires labels we don't produce yet -- namely, those at the boundary between PDcoord~0.9-1 and SRLM, and between PDcoord~0.9-1 and BG +} + desc_io = { "hipp": "equivol" if "equivolume" in config["laminar_coords_method"] else "laplace", "dentate": "laplace", @@ -25,109 +42,90 @@ ruleorder: atlas_label_to_unfold_nii > atlas_metric_to_unfold_nii #temporary un # --- isosurface generation --- -rule prep_hipp_coords_for_meshing: +rule get_label_mask: input: - coords=lambda wildcards: bids( - root=work, - datatype="coords", - dir="IO", - label="{label}", - suffix="coords.nii.gz", - desc=desc_io[wildcards.label], - space="corobl", - hemi="{hemi}", - **inputs.subj_wildcards - ), labelmap=get_labels_for_laplace, params: - gm_labels=lambda wildcards: config["laplace_labels"]["IO"]["gm"], - gm_noDG_labels=lambda wildcards: config["laplace_labels"]["IO"]["gm_noDG"], - src_labels=lambda wildcards: config["laplace_labels"]["IO"]["src"], - sink_labels=lambda wildcards: config["laplace_labels"]["IO"]["sink"], - output: - coords=temp( - bids( - root=work, - datatype="surf", - suffix="coords.nii.gz", - space="corobl", - desc="formesh", - hemi="{hemi}", - label="{label,hipp}", - **inputs.subj_wildcards - ) + gm_labels=lambda wildcards: " ".join( + [str(lbl) for lbl in gm_labels[wildcards.label]] ), + output: mask=temp( bids( root=work, - datatype="surf", + datatype="anat", suffix="mask.nii.gz", space="corobl", - desc="formesh", + desc="GM", hemi="{hemi}", - label="{label,hipp}", + label="{label}", **inputs.subj_wildcards ) ), container: config["singularity"]["autotop"] - script: - "../scripts/prep_hipp_coords_for_meshing.py" + shell: + "c3d -background -1 {input} -retain-labels {params} -binarize {output}" -rule prep_dentate_coords_for_meshing: +rule get_sink_mask: input: - coords=lambda wildcards: bids( - root=work, - datatype="coords", - dir="IO", - label="{label}", - suffix="coords.nii.gz", - desc=desc_io[wildcards.label], - space="corobl", - hemi="{hemi}", - **inputs.subj_wildcards - ), labelmap=get_labels_for_laplace, params: - gm_labels=lambda wildcards: config["laplace_labels"]["IO"]["gm"], + labels=lambda wildcards: " ".join( + [str(lbl) for lbl in sink_labels[wildcards.label]] + ), output: - coords=temp( + mask=temp( bids( root=work, - datatype="surf", - suffix="coords.nii.gz", + datatype="anat", + suffix="mask.nii.gz", space="corobl", - desc="formesh", + desc="sink", hemi="{hemi}", - label="{label,dentate}", + label="{label}", **inputs.subj_wildcards ) ), + container: + config["singularity"]["autotop"] + shell: + "c3d {input} -background -1 -retain-labels {params} -binarize {output}" + + +rule get_src_mask: + input: + labelmap=get_labels_for_laplace, + params: + labels=lambda wildcards: " ".join( + [str(lbl) for lbl in src_labels[wildcards.label]] + ), + output: mask=temp( bids( root=work, - datatype="surf", + datatype="anat", suffix="mask.nii.gz", space="corobl", - desc="formesh", + desc="src", hemi="{hemi}", - label="{label,dentate}", + label="{label}", **inputs.subj_wildcards ) ), container: config["singularity"]["autotop"] - script: - "../scripts/prep_dentate_coords_for_meshing.py" + shell: + "c3d {input} -background -1 -retain-labels {params} -binarize {output}" -rule get_label_mask: +rule get_nan_mask: input: labelmap=get_labels_for_laplace, params: - gm_labels=lambda wildcards: " ".join( - [str(lbl) for lbl in gm_labels[wildcards.label]] + labels=lambda wildcards: " ".join( + [str(lbl) for lbl in nan_labels[wildcards.label]] ), output: mask=temp( @@ -136,7 +134,7 @@ rule get_label_mask: datatype="anat", suffix="mask.nii.gz", space="corobl", - desc="GM", + desc="nan", hemi="{hemi}", label="{label}", **inputs.subj_wildcards @@ -145,12 +143,12 @@ rule get_label_mask: container: config["singularity"]["autotop"] shell: - "c3d {input} -retain-labels {params} -binarize {output}" + "c3d {input} -background -1 -retain-labels {params} -binarize {output}" rule gen_native_mesh: input: - nii=lambda wildcards: bids( + coords=lambda wildcards: bids( root=work, datatype="coords", dir="IO", @@ -161,12 +159,32 @@ rule gen_native_mesh: hemi="{hemi}", **inputs.subj_wildcards ), - mask=bids( + nan_mask=bids( root=work, datatype="anat", suffix="mask.nii.gz", space="corobl", - desc="GM", + desc="nan", + hemi="{hemi}", + label="{label}", + **inputs.subj_wildcards + ), + sink_mask=bids( + root=work, + datatype="anat", + suffix="mask.nii.gz", + space="corobl", + desc="sink", + hemi="{hemi}", + label="{label}", + **inputs.subj_wildcards + ), + src_mask=bids( + root=work, + datatype="anat", + suffix="mask.nii.gz", + space="corobl", + desc="src", hemi="{hemi}", label="{label}", **inputs.subj_wildcards @@ -283,14 +301,15 @@ rule warp_native_mesh_to_unfold: rule compute_halfthick_mask: input: - coords=bids( + coords=lambda wildcards: bids( root=work, - datatype="surf", + datatype="coords", + dir="IO", + label="{label}", suffix="coords.nii.gz", + desc=desc_io[wildcards.label], space="corobl", - desc="formesh", hemi="{hemi}", - label="{label}", **inputs.subj_wildcards ), mask=bids( diff --git a/hippunfold/workflow/scripts/gen_isosurface.py b/hippunfold/workflow/scripts/gen_isosurface.py index ecf66074..0ee9ff6c 100644 --- a/hippunfold/workflow/scripts/gen_isosurface.py +++ b/hippunfold/workflow/scripts/gen_isosurface.py @@ -54,29 +54,73 @@ def remove_nan_points_faces(vertices, faces): return (new_vertices, new_faces) -# Load the NIfTI file -nifti_img = nib.load(snakemake.input.nii) -data = nifti_img.get_fdata() +from scipy.ndimage import binary_dilation -# Load the mask -mask_img = nib.load(snakemake.input.mask) -mask = mask_img.get_fdata() +def get_adjacent_voxels(mask_a, mask_b): + """ + Create a mask for voxels where label A is adjacent to label B. + Parameters: + - mask_a (np.ndarray): A 3D binary mask for label A. + - mask_b (np.ndarray): A 3D binary mask for label B. -affine = nifti_img.affine + Returns: + - np.ndarray: A 3D mask where adjacent voxels for label A and label B are marked as True. + """ + # Dilate each mask to identify neighboring regions + dilated_a = binary_dilation(mask_a) + dilated_b = binary_dilation(mask_b) + + # Find adjacency: voxels of A touching B and B touching A + adjacency_mask = (dilated_a.astype("bool") & mask_b.astype("bool")) | ( + dilated_b.astype("bool") & mask_a.astype("bool") + ) + + return adjacency_mask + + +# Load the coords image +coords_img = nib.load(snakemake.input.coords) +coords = coords_img.get_fdata() + + +# Load the nan mask +nan_mask_img = nib.load(snakemake.input.nan_mask) +nan_mask = nan_mask_img.get_fdata() + +# Load the sink mask +sink_mask_img = nib.load(snakemake.input.sink_mask) +sink_mask = sink_mask_img.get_fdata() + +# Load the src mask +src_mask_img = nib.load(snakemake.input.src_mask) +src_mask = src_mask_img.get_fdata() + + +affine = coords_img.affine # Get voxel spacing from the header -voxel_spacing = nifti_img.header.get_zooms()[:3] +voxel_spacing = coords_img.header.get_zooms()[:3] # Create a PyVista grid grid = pv.ImageData() -grid.dimensions = np.array(data.shape) + 1 # Add 1 to include the boundary voxels +grid.dimensions = np.array(coords.shape) + 1 # Add 1 to include the boundary voxels grid.spacing = (1, 1, 1) # Use unit spacing and zero origin since we will apply affine grid.origin = (0, 0, 0) +# update the coords data to add the nans and sink +coords[nan_mask == 1] = np.nan +coords[sink_mask == 1] = 1.1 # since sink being zero creates a false boundary + +# we also need to use a nan mask for the voxels where src and sink meet directly +# (since this is another false boundary).. +src_sink_nan_mask = get_adjacent_voxels(sink_mask, src_mask) +coords[src_sink_nan_mask == 1] = np.nan + + # Add the scalar field -grid.cell_data["values"] = np.where(mask == 0, np.nan, data).flatten(order="F") +grid.cell_data["values"] = coords.flatten(order="F") grid = grid.cells_to_points("values") # apply the affine diff --git a/hippunfold/workflow/scripts/prep_dentate_coords_for_meshing.py b/hippunfold/workflow/scripts/prep_dentate_coords_for_meshing.py deleted file mode 100644 index c47e56b7..00000000 --- a/hippunfold/workflow/scripts/prep_dentate_coords_for_meshing.py +++ /dev/null @@ -1,39 +0,0 @@ -import nibabel as nib -import numpy as np -from scipy.ndimage import binary_dilation -from astropy.convolution import convolve - -coords_nib = nib.load(snakemake.input.coords) -label_nib = nib.load(snakemake.input.labelmap) - -coords = coords_nib.get_fdata() -label = label_nib.get_fdata() - -gm_mask = np.zeros(label.shape, dtype=bool) -for lbl in snakemake.params.gm_labels: - gm_mask = np.where(label == lbl, True, gm_mask) - - -# set everywhere outside coords mask to be nan, so they don't contribute to extrapolation -coords_withnan = np.where(gm_mask, coords, np.nan) - -# perform convolution with astropy to fill in nans -hl = np.ones([3, 3, 3]) -hl = hl / np.sum(hl) -convolved_coords = convolve( - coords_withnan, hl, nan_treatment="interpolate", preserve_nan=False -) - -# dilating the mask into the src and sink -structuring_element = np.ones((3, 3, 3), dtype=bool) -dilated_mask = binary_dilation(gm_mask, structuring_element) - -# for dentate, we use the convolved coords (hopefully the smoothing is helpful here for making a surface - -# save the extrapolated coords and the updated mask -nib.Nifti1Image(convolved_coords, coords_nib.affine, coords_nib.header).to_filename( - snakemake.output.coords -) -nib.Nifti1Image(dilated_mask, coords_nib.affine, coords_nib.header).to_filename( - snakemake.output.mask -) diff --git a/hippunfold/workflow/scripts/prep_hipp_coords_for_meshing.py b/hippunfold/workflow/scripts/prep_hipp_coords_for_meshing.py deleted file mode 100644 index 4ecc38d8..00000000 --- a/hippunfold/workflow/scripts/prep_hipp_coords_for_meshing.py +++ /dev/null @@ -1,57 +0,0 @@ -import nibabel as nib -import numpy as np -from scipy.ndimage import binary_dilation -from astropy.convolution import convolve - -coords_nib = nib.load(snakemake.input.coords) -label_nib = nib.load(snakemake.input.labelmap) - -coords = coords_nib.get_fdata() -label = label_nib.get_fdata() - -# get the source and sink regions -src_mask = np.zeros(label.shape, dtype=bool) -for lbl in snakemake.params.src_labels: - src_mask = np.where(label == lbl, True, src_mask) - -sink_mask = np.zeros(label.shape, dtype=bool) -for lbl in snakemake.params.sink_labels: - sink_mask = np.where(label == lbl, True, sink_mask) - -gm_mask = np.zeros(label.shape, dtype=bool) -for lbl in snakemake.params.gm_labels: - gm_mask = np.where(label == lbl, True, gm_mask) - -gm_noDG_mask = np.zeros(label.shape, dtype=bool) -for lbl in snakemake.params.gm_noDG_labels: - gm_noDG_mask = np.where(label == lbl, True, gm_noDG_mask) - - -# set everywhere outside coords mask to be nan, so they don't contribute to extrapolation -coords_withnan = np.where(gm_mask, coords, np.nan) - -# perform convolution with astropy to fill in nans -hl = np.ones([3, 3, 3]) -hl = hl / np.sum(hl) -convolved_coords = convolve( - coords_withnan, hl, nan_treatment="interpolate", preserve_nan=False -) - -# dilating the mask into the src and sink -structuring_element = np.ones((3, 3, 3), dtype=bool) -dilated_mask = binary_dilation(gm_noDG_mask, structuring_element) & ( - src_mask | sink_mask -) -updated_mask = np.where(dilated_mask, 1, gm_noDG_mask) - -# get the updated coords: -updated_coords = np.where(gm_mask, coords, np.where(dilated_mask, convolved_coords, 0)) - - -# save the extrapolated coords and the updated mask -nib.Nifti1Image(convolved_coords, coords_nib.affine, coords_nib.header).to_filename( - snakemake.output.coords -) -nib.Nifti1Image(updated_mask, coords_nib.affine, coords_nib.header).to_filename( - snakemake.output.mask -)