orthogonalize AP/PD src/sink edges

hippunfold/workflow/rules/download.smk

@@ -120,6 +120,8 @@ rule import_template_dseg_dentate:
         "subj"
     container:
         config["singularity"]["autotop"]
+    conda:
+        "../envs/c3d.yaml"
     shell:
         "{params.copy_or_flip_cmd} {output.template_seg}"
 

hippunfold/workflow/rules/native_surf.smk

@@ -188,7 +188,7 @@ rule get_boundary_vertices:
     output:
         label_gii=bids(
             root=work,
-            datatype="surf",
+            datatype="coords",
             suffix="boundary.label.gii",
             space="corobl",
             hemi="{hemi}",
@@ -231,7 +231,7 @@ rule map_src_sink_sdt_to_surf:
     output:
         sdt=bids(
             root=work,
-            datatype="surf",
+            datatype="coords",
             suffix="sdt.shape.gii",
             space="corobl",
             hemi="{hemi}",
@@ -250,54 +250,129 @@ rule map_src_sink_sdt_to_surf:
         "wb_command -volume-to-surface-mapping {input.sdt} {input.surf_gii} {output.sdt} -trilinear"
 
 
-rule laplace_beltrami:
+rule postproc_boundary_vertices:
+    """ ensures non-overlapping and full labelling of AP/PD edges """
     input:
-        surf_gii=bids(
-            root=root,
-            datatype="surf",
-            suffix="midthickness.surf.gii",
+        ap_src=bids(
+            root=work,
+            datatype="coords",
+            suffix="sdt.shape.gii",
             space="corobl",
             hemi="{hemi}",
+            dir="AP",
+            desc="src",
             label="{label}",
             **inputs.subj_wildcards,
         ),
-        src_sdt=bids(
+        ap_sink=bids(
             root=work,
-            datatype="surf",
+            datatype="coords",
             suffix="sdt.shape.gii",
             space="corobl",
             hemi="{hemi}",
-            dir="{dir}",
+            dir="AP",
+            desc="sink",
+            label="{label}",
+            **inputs.subj_wildcards,
+        ),
+        pd_src=bids(
+            root=work,
+            datatype="coords",
+            suffix="sdt.shape.gii",
+            space="corobl",
+            hemi="{hemi}",
+            dir="PD",
             desc="src",
             label="{label}",
             **inputs.subj_wildcards,
         ),
-        sink_sdt=bids(
+        pd_sink=bids(
             root=work,
-            datatype="surf",
+            datatype="coords",
             suffix="sdt.shape.gii",
             space="corobl",
             hemi="{hemi}",
-            dir="{dir}",
+            dir="PD",
             desc="sink",
             label="{label}",
             **inputs.subj_wildcards,
         ),
-        boundary=bids(
+        edges=bids(
             root=work,
-            datatype="surf",
+            datatype="coords",
             suffix="boundary.label.gii",
             space="corobl",
             hemi="{hemi}",
             label="{label}",
             **inputs.subj_wildcards,
         ),
     params:
-        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"
+        min_terminal_vertices=5,  # min number of vertices per src/sink
+        max_iterations=100,
+        shifting_epsilon=0.1,  #could be proportional to voxel spacing
+    output:
+        ap=bids(
+            root=work,
+            datatype="coords",
+            suffix="mask.label.gii",
+            space="corobl",
+            hemi="{hemi}",
+            dir="AP",
+            desc="srcsink",
+            label="{label}",
+            **inputs.subj_wildcards,
+        ),
+        pd=bids(
+            root=work,
+            datatype="coords",
+            suffix="mask.label.gii",
+            space="corobl",
+            hemi="{hemi}",
+            dir="PD",
+            desc="srcsink",
+            label="{label}",
+            **inputs.subj_wildcards,
+        ),
+    container:
+        config["singularity"]["autotop"]
+    log:
+        bids(
+            root="logs",
+            datatype="postproc_boundary_vertices",
+            suffix="log.txt",
+            hemi="{hemi}",
+            label="{label}",
+            **inputs.subj_wildcards,
+        ),
+    conda:
+        "../envs/pyvista.yaml"
+    group:
+        "subj"
+    script:
+        "../scripts/postproc_boundary_vertices.py"
+
+
+rule laplace_beltrami:
+    input:
+        surf_gii=bids(
+            root=root,
+            datatype="surf",
+            suffix="midthickness.surf.gii",
+            space="corobl",
+            hemi="{hemi}",
+            label="{label}",
+            **inputs.subj_wildcards,
+        ),
+        src_sink_mask=bids(
+            root=work,
+            datatype="coords",
+            suffix="mask.label.gii",
+            space="corobl",
+            hemi="{hemi}",
+            dir="{dir}",
+            desc="srcsink",
+            label="{label}",
+            **inputs.subj_wildcards,
         ),
     output:
         coords=bids(
@@ -311,6 +386,16 @@ rule laplace_beltrami:
             label="{label}",
             **inputs.subj_wildcards,
         ),
+    log:
+        bids(
+            root="logs",
+            datatype="laplace_beltrami",
+            suffix="log.txt",
+            hemi="{hemi}",
+            dir="{dir}",
+            label="{label}",
+            **inputs.subj_wildcards,
+        ),
     group:
         "subj"
     container:

hippunfold/workflow/scripts/get_boundary_vertices.py

@@ -2,6 +2,7 @@
 import numpy as np
 import nibabel as nib
 import nibabel.gifti as gifti
+from nibabel.gifti.gifti import intent_codes
 from collections import defaultdict
 from lib.utils import setup_logger
 
@@ -50,11 +51,24 @@ def read_surface_from_gifti(surf_gii):
     faces = surf.agg_data("NIFTI_INTENT_TRIANGLE")
     faces_pv = np.hstack([np.full((faces.shape[0], 1), 3), faces])  # PyVista format
 
-    return pv.PolyData(vertices, faces_pv)
+    # Find the first darray that represents vertices (NIFTI_INTENT_POINTSET)
+    vertices_darray = next(
+        (
+            darray
+            for darray in surf.darrays
+            if darray.intent == intent_codes["NIFTI_INTENT_POINTSET"]
+        ),
+        None,
+    )
+
+    # Extract metadata as a dictionary (return empty dict if no metadata)
+    metadata = dict(vertices_darray.meta) if vertices_darray else {}
+
+    return pv.PolyData(vertices, faces_pv), metadata
 
 
 logger.info("Loading surface from GIFTI...")
-surface = read_surface_from_gifti(snakemake.input.surf_gii)
+surface, metadata = read_surface_from_gifti(snakemake.input.surf_gii)
 logger.info(f"Surface loaded: {surface.n_points} vertices, {surface.n_faces} faces.")
 
 
@@ -67,6 +81,33 @@ def read_surface_from_gifti(surf_gii):
     f"Boundary scalar array created. {np.sum(boundary_scalars)} boundary vertices marked."
 )
 
+
+# Find the largest connected component within this sub-mesh
+logger.info("Applying largest connected components")
+
+# Extract points that are within the boundary scalars
+sub_mesh = pv.PolyData(surface.points, surface.faces).extract_points(
+    boundary_scalars.astype(bool), adjacent_cells=True
+)
+
+# Compute connectivity to find the largest connected component
+connected_sub_mesh = sub_mesh.connectivity("largest")
+
+# Get indices of the largest component in the sub-mesh
+largest_component_mask = (
+    connected_sub_mesh.point_data["RegionId"] == 0
+)  # Largest component has RegionId 0
+largest_component_indices = connected_sub_mesh.point_data["vtkOriginalPointIds"][
+    largest_component_mask
+]
+
+# Create an array for all points in the original surface
+boundary_scalars = np.zeros(surface.n_points, dtype=np.int32)
+
+# Keep only the largest component
+boundary_scalars[largest_component_indices] = 1
+
+
 logger.info("Saving GIFTI label file...")
 
 # Create a GIFTI label data array
@@ -92,6 +133,10 @@ def read_surface_from_gifti(surf_gii):
 # Assign label table to GIFTI image
 gii_img = gifti.GiftiImage(darrays=[gii_data], labeltable=label_table)
 
+# set structure metadata
+gii_img.meta["AnatomicalStructurePrimary"] = metadata["AnatomicalStructurePrimary"]
+
+
 # Save the label file
 gii_img.to_filename(snakemake.output.label_gii)
 logger.info(f"GIFTI label file saved as '{snakemake.output.label_gii}'.")