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

@@ -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(

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