minor bug fixes, mostly formatting

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,12 +250,12 @@ rule map_src_sink_sdt_to_surf:
         "wb_command -volume-to-surface-mapping {input.sdt} {input.surf_gii} {output.sdt} -trilinear"
 
 
-rule joint_smooth_ap_pd_edges:
+rule postproc_boundary_vertices:
     """ ensures non-overlapping and full labelling of AP/PD edges """
     input:
         ap_src=bids(
             root=work,
-            datatype="surf",
+            datatype="coords",
             suffix="sdt.shape.gii",
             space="corobl",
             hemi="{hemi}",
@@ -266,7 +266,7 @@ rule joint_smooth_ap_pd_edges:
         ),
         ap_sink=bids(
             root=work,
-            datatype="surf",
+            datatype="coords",
             suffix="sdt.shape.gii",
             space="corobl",
             hemi="{hemi}",
@@ -277,7 +277,7 @@ rule joint_smooth_ap_pd_edges:
         ),
         pd_src=bids(
             root=work,
-            datatype="surf",
+            datatype="coords",
             suffix="sdt.shape.gii",
             space="corobl",
             hemi="{hemi}",
@@ -288,7 +288,7 @@ rule joint_smooth_ap_pd_edges:
         ),
         pd_sink=bids(
             root=work,
-            datatype="surf",
+            datatype="coords",
             suffix="sdt.shape.gii",
             space="corobl",
             hemi="{hemi}",
@@ -299,7 +299,7 @@ rule joint_smooth_ap_pd_edges:
         ),
         edges=bids(
             root=work,
-            datatype="surf",
+            datatype="coords",
             suffix="boundary.label.gii",
             space="corobl",
             hemi="{hemi}",
@@ -311,7 +311,7 @@ rule joint_smooth_ap_pd_edges:
     output:
         ap=bids(
             root=work,
-            datatype="surf",
+            datatype="coords",
             suffix="mask.label.gii",
             space="corobl",
             hemi="{hemi}",
@@ -322,7 +322,7 @@ rule joint_smooth_ap_pd_edges:
         ),
         pd=bids(
             root=work,
-            datatype="surf",
+            datatype="coords",
             suffix="mask.label.gii",
             space="corobl",
             hemi="{hemi}",
@@ -354,7 +354,7 @@ rule laplace_beltrami:
         ),
         src_sink_mask=bids(
             root=work,
-            datatype="surf",
+            datatype="coords",
             suffix="mask.label.gii",
             space="corobl",
             hemi="{hemi}",

hippunfold/workflow/scripts/get_boundary_vertices.py

@@ -67,14 +67,31 @@ 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
 )
-largest_component_indices = sub_mesh.connectivity(largest=True).point_data["RegionId"]
-boundary_scalars = np.zeros(surface.n_points, dtype=int)
-boundary_scalars[largest_component.point_data["vtkOriginalPointIds"]] = 1
-logger.info("Applying largest connected components")
+
+# 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...")

hippunfold/workflow/scripts/laplace_beltrami.py

@@ -105,7 +105,9 @@ def solve_laplace_beltrami_open_mesh(vertices, faces, boundary_conditions=None):
 src_vals = [0 for i in range(len(src_indices))]
 sink_vals = [1 for i in range(len(sink_indices))]
 
-boundary_conditions = dict(zip(src_indices + sink_indices, src_vals + sink_vals))
+boundary_conditions = dict(
+    zip(list(src_indices) + list(sink_indices), src_vals + sink_vals)
+)
 
 
 coords = solve_laplace_beltrami_open_mesh(vertices, faces, boundary_conditions)

hippunfold/workflow/scripts/postproc_boundary_vertices.py

@@ -8,20 +8,22 @@
 nmin = snakemake.params.min_terminal_vertices
 
 logger.info("Loading surface from GIFTI...")
+
 edges = nib.load(snakemake.input.edges).agg_data()
 ap_src = nib.load(snakemake.input.ap_src).agg_data()
 ap_sink = nib.load(snakemake.input.ap_sink).agg_data()
 pd_src = nib.load(snakemake.input.pd_src).agg_data()
 pd_sink = nib.load(snakemake.input.pd_sink).agg_data()
 
-distances = np.vstack(
-    (ap_src[edges == 1], ap_sink[edges == 1], pd_src[edges == 1], pd_sink[edges == 1])
-)
-num_vertices = distances.shape[0]
-
 
 logger.info("Assigning labels apsrc, apsink, pdsrc, pdsink")
 
+distances = np.vstack(
+    (ap_src[edges == 1], ap_sink[edges == 1], pd_src[edges == 1], pd_sink[edges == 1])
+).T
+scaling_factors = np.ones((4))
+num_labels = 4  # starts at 0
+
 max_iterations = 10  # Prevent infinite loops
 for _ in range(max_iterations):
     # Scale distances
@@ -43,19 +45,24 @@
     # Update scaling factors for underrepresented labels
     for k in range(num_labels):
         if label_counts[k] < nmin:
-            scaling_factors[k] *= 1.5  # Increase competitiveness of the label
+            scaling_factors[k] *= 0.9  # Increase competitiveness of the label
 
 # Ensure all labels are represented
-logger.info("Final label counts:", label_counts)
-
-ap_srcsink = np.zeros((len(edges)))
-ap_srcsink[edges == 1][labels == 0] = 1
-ap_srcsink[edges == 1][labels == 1] = 2
-gii_img = gifti.GiftiImage(darrays=[ap_srcsink])
-nib.save(gii_img, snakemake.outputs.ap)
-
-pd_srcsink = np.zeros((len(edges)))
-pd_srcsink[edges == 1][labels == 2] = 1
-pd_srcsink[edges == 1][labels == 3] = 2
-gii_img = gifti.GiftiImage(darrays=[pd_srcsink])
-nib.save(gii_img, snakemake.outputs.pd)
+logger.info(["Final label counts:", label_counts])
+
+ap_srcsink = np.zeros((len(edges)), dtype=np.int32)
+idx_edges = np.where(edges == 1)[0]
+ap_srcsink[idx_edges[labels == 0]] = 1
+ap_srcsink[idx_edges[labels == 1]] = 2
+gii_img = gifti.GiftiImage(
+    darrays=[gifti.GiftiDataArray(ap_srcsink, intent="NIFTI_INTENT_LABEL")]
+)
+nib.save(gii_img, snakemake.output.ap)
+
+pd_srcsink = np.zeros((len(edges)), dtype=np.int32)
+pd_srcsink[idx_edges[labels == 2]] = 1
+pd_srcsink[idx_edges[labels == 3]] = 2
+gii_img = gifti.GiftiImage(
+    darrays=[gifti.GiftiDataArray(pd_srcsink, intent="NIFTI_INTENT_LABEL")]
+)
+nib.save(gii_img, snakemake.output.pd)