From 64546f9a4d2a7d726e3ea6dc54abb667ce55cf3a Mon Sep 17 00:00:00 2001 From: Ali Khan Date: Mon, 10 Feb 2025 16:52:50 -0500 Subject: [PATCH 1/5] refactoring of the laplace beltrami rule (#371) - splits into different rules for extracting the boundary, and getting distances to src/sink masks - those new rules now produce surface maps (for easier introspection), which are passed along to the laplace beltrami rule - the dentate AP src in bbhist hemi-L was problematic with existing approach, since the PD src/sink were always closer than the AP src, so no vertices were being labelled with AP src. - in this version, we do things a little different: 1. AP and PD are done independently now 2. Instead of nearest neighbor, we use a distance threshold, which is defined based on a signed distance transform of the mask, mapped to the surface. 3. Actually there is a minimum distance, and maximum distance, along with the minimum number of vertices required in the src/sink. The threshold is swept until the desired number is reached, to deal with cases such as the above. - TODO: the get_boundary_vertices script could be further improved to perform connected components and picking the largest one, to avoid being affected by holes or other defects.. - TODO: the laplace beltrami solver doesn't take as long to run now that I am decimating alot, but if we want to further optimize things, I think the step of altering the laplacian based on the boundary conditions, after it has been made a sparse matrix, is where the inefficiency lies (based on profiling). Could try to avoid setting those weights in the first place instead of setting it after the fact.. - NOTE: solve_laplace_beltrami_open_mesh() is kept as is (except for some logging statements) --- hippunfold/workflow/rules/coords.smk | 45 +++--- hippunfold/workflow/rules/native_surf.smk | 119 ++++++++++++++-- .../workflow/scripts/get_boundary_vertices.py | 97 +++++++++++++ .../workflow/scripts/laplace_beltrami.py | 132 ++++++++---------- 4 files changed, 285 insertions(+), 108 deletions(-) create mode 100644 hippunfold/workflow/scripts/get_boundary_vertices.py diff --git a/hippunfold/workflow/rules/coords.smk b/hippunfold/workflow/rules/coords.smk index dfcd2529..8db4cb09 100644 --- a/hippunfold/workflow/rules/coords.smk +++ b/hippunfold/workflow/rules/coords.smk @@ -21,21 +21,13 @@ def get_gm_labels(wildcards): return lbl_list -def get_sink_labels(wildcards): +def get_src_sink_labels(wildcards): lbl_list = " ".join( [ str(lbl) - for lbl in config["laplace_labels"][wildcards.label][wildcards.dir]["sink"] - ] - ) - return lbl_list - - -def get_src_labels(wildcards): - lbl_list = " ".join( - [ - str(lbl) - for lbl in config["laplace_labels"][wildcards.label][wildcards.dir]["src"] + for lbl in config["laplace_labels"][wildcards.label][wildcards.dir][ + wildcards.srcsink + ] ] ) return lbl_list @@ -100,11 +92,11 @@ def get_inputs_laplace(wildcards): return files -rule get_sink_mask: +rule get_src_sink_mask: input: labelmap=get_labels_for_laplace, params: - labels=get_sink_labels, + labels=get_src_sink_labels, output: mask=bids( root=work, @@ -112,7 +104,7 @@ rule get_sink_mask: suffix="mask.nii.gz", space="corobl", dir="{dir}", - desc="sink", + desc="{srcsink,src|sink}", hemi="{hemi}", label="{label}", **inputs.subj_wildcards, @@ -127,19 +119,28 @@ rule get_sink_mask: "c3d {input} -background -1 -retain-labels {params} -binarize {output}" -rule get_src_mask: +rule get_src_sink_sdt: + """calculate signed distance transform (negative inside, positive outside)""" input: - labelmap=get_labels_for_laplace, - params: - labels=get_src_labels, - output: mask=bids( root=work, datatype="coords", suffix="mask.nii.gz", space="corobl", dir="{dir}", - desc="src", + desc="{srcsink}", + hemi="{hemi}", + label="{label}", + **inputs.subj_wildcards, + ), + output: + sdt=bids( + root=work, + datatype="coords", + suffix="sdt.nii.gz", + space="corobl", + dir="{dir}", + desc="{srcsink}", hemi="{hemi}", label="{label}", **inputs.subj_wildcards, @@ -151,7 +152,7 @@ rule get_src_mask: group: "subj" shell: - "c3d {input} -background -1 -retain-labels {params} -binarize {output}" + "c3d {input} -sdt -o {output}" rule get_nan_mask: diff --git a/hippunfold/workflow/rules/native_surf.smk b/hippunfold/workflow/rules/native_surf.smk index 38d87633..3b1f1dda 100644 --- a/hippunfold/workflow/rules/native_surf.smk +++ b/hippunfold/workflow/rules/native_surf.smk @@ -66,7 +66,7 @@ rule gen_native_mesh: params: threshold=lambda wildcards: surf_thresholds[wildcards.surfname], decimate_opts={ - "reduction": 0.5, + "reduction": 0.9, "feature_angle": 25, "preserve_topology": True, }, @@ -173,7 +173,7 @@ rule smooth_surface: # --- creating unfold surface from native anatomical, including post-processing -rule laplace_beltrami: +rule get_boundary_vertices: input: surf_gii=bids( root=root, @@ -184,25 +184,122 @@ rule laplace_beltrami: label="{label}", **inputs.subj_wildcards, ), - seg=get_labels_for_laplace, - params: - srcsink_labels=lambda wildcards: config["laplace_labels"][wildcards.label], output: - coords_AP=bids( + label_gii=bids( + root=work, + datatype="surf", + suffix="boundary.label.gii", + space="corobl", + hemi="{hemi}", + label="{label}", + **inputs.subj_wildcards, + ), + group: + "subj" + container: + config["singularity"]["autotop"] + conda: + "../envs/pyvista.yaml" + script: + "../scripts/get_boundary_vertices.py" + + +rule map_src_sink_sdt_to_surf: + """ Maps the distance to src/sink mask """ + input: + surf_gii=bids( + root=root, + datatype="surf", + suffix="midthickness.surf.gii", + space="corobl", + hemi="{hemi}", + label="{label}", + **inputs.subj_wildcards, + ), + sdt=bids( root=work, datatype="coords", - dir="AP", - suffix="coords.shape.gii", - desc="laplace", + suffix="sdt.nii.gz", space="corobl", + dir="{dir}", + desc="{srcsink}", hemi="{hemi}", label="{label}", **inputs.subj_wildcards, ), - coords_PD=bids( + output: + sdt=bids( + root=work, + datatype="surf", + suffix="sdt.shape.gii", + space="corobl", + hemi="{hemi}", + dir="{dir}", + desc="{srcsink}", + label="{label}", + **inputs.subj_wildcards, + ), + container: + config["singularity"]["autotop"] + conda: + "../envs/workbench.yaml" + group: + "subj" + shell: + "wb_command -volume-to-surface-mapping {input.sdt} {input.surf_gii} {output.sdt} -trilinear" + + +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_sdt=bids( + root=work, + datatype="surf", + suffix="sdt.shape.gii", + space="corobl", + hemi="{hemi}", + dir="{dir}", + desc="src", + label="{label}", + **inputs.subj_wildcards, + ), + sink_sdt=bids( + root=work, + datatype="surf", + suffix="sdt.shape.gii", + space="corobl", + hemi="{hemi}", + dir="{dir}", + desc="sink", + label="{label}", + **inputs.subj_wildcards, + ), + boundary=bids( + root=work, + datatype="surf", + suffix="boundary.label.gii", + space="corobl", + hemi="{hemi}", + label="{label}", + **inputs.subj_wildcards, + ), + params: + min_dist_threshold=0.3, + max_dist_threshold=1, + min_terminal_vertices=5, + output: + coords=bids( root=work, datatype="coords", - dir="PD", + dir="{dir}", suffix="coords.shape.gii", desc="laplace", space="corobl", diff --git a/hippunfold/workflow/scripts/get_boundary_vertices.py b/hippunfold/workflow/scripts/get_boundary_vertices.py new file mode 100644 index 00000000..1e5101d3 --- /dev/null +++ b/hippunfold/workflow/scripts/get_boundary_vertices.py @@ -0,0 +1,97 @@ +import pyvista as pv +import numpy as np +import nibabel as nib +import nibabel.gifti as gifti +from collections import defaultdict +from lib.utils import setup_logger + +# Setup logger +log_file = snakemake.log[0] if snakemake.log else None +logger = setup_logger(log_file) + + +def find_boundary_vertices(mesh): + """ + Find boundary vertices of a 3D mesh. + + Args: + mesh + Returns: + list: List of vertex indices that are boundary vertices, sorted in ascending order. + """ + vertices = mesh.points + faces = mesh.faces.reshape((-1, 4))[:, 1:4] # Extract triangle indices + + edge_count = defaultdict(int) + # Step 1: Count edge occurrences + for face in faces: + # Extract edges from the face, ensure consistent ordering (min, max) + edges = [ + tuple(sorted((face[0], face[1]))), + tuple(sorted((face[1], face[2]))), + tuple(sorted((face[2], face[0]))), + ] + for edge in edges: + edge_count[edge] += 1 + # Step 2: Identify boundary edges + boundary_edges = [edge for edge, count in edge_count.items() if count == 1] + # Step 3: Collect boundary vertices + boundary_vertices = set() + for edge in boundary_edges: + boundary_vertices.update(edge) + # Convert the set to a sorted list (array) + return np.array(sorted(boundary_vertices), dtype=np.int32) + + +def read_surface_from_gifti(surf_gii): + """Load a surface mesh from a GIFTI file.""" + surf = nib.load(surf_gii) + vertices = surf.agg_data("NIFTI_INTENT_POINTSET") + 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) + + +logger.info("Loading surface from GIFTI...") +surface = read_surface_from_gifti(snakemake.input.surf_gii) +logger.info(f"Surface loaded: {surface.n_points} vertices, {surface.n_faces} faces.") + + +logger.info("Find boundary vertices") +boundary_indices = find_boundary_vertices(surface) + +boundary_scalars = np.zeros(surface.n_points, dtype=np.int32) # Default is 0 +boundary_scalars[boundary_indices] = 1 # Set boundary vertices to 1 +logger.info( + f"Boundary scalar array created. {np.sum(boundary_scalars)} boundary vertices marked." +) + +logger.info("Saving GIFTI label file...") + +# Create a GIFTI label data array +gii_data = gifti.GiftiDataArray(boundary_scalars, intent="NIFTI_INTENT_LABEL") + +# Create a Label Table (LUT) +label_table = gifti.GiftiLabelTable() + +# Define Background label (key 0) +background_label = gifti.GiftiLabel( + key=0, red=1.0, green=1.0, blue=1.0, alpha=0.0 +) # Transparent +background_label.label = "Background" +label_table.labels.append(background_label) + +# Define Boundary label (key 1) +boundary_label = gifti.GiftiLabel( + key=1, red=1.0, green=0.0, blue=0.0, alpha=1.0 +) # Red color +boundary_label.label = "Boundary" +label_table.labels.append(boundary_label) + +# Assign label table to GIFTI image +gii_img = gifti.GiftiImage(darrays=[gii_data], labeltable=label_table) + +# Save the label file +gii_img.to_filename(snakemake.output.label_gii) +logger.info(f"GIFTI label file saved as '{snakemake.output.label_gii}'.") diff --git a/hippunfold/workflow/scripts/laplace_beltrami.py b/hippunfold/workflow/scripts/laplace_beltrami.py index de687ff1..e9cfabc8 100644 --- a/hippunfold/workflow/scripts/laplace_beltrami.py +++ b/hippunfold/workflow/scripts/laplace_beltrami.py @@ -1,41 +1,26 @@ import numpy as np -from scipy.sparse import diags, linalg, lil_matrix import nibabel as nib -from scipy.interpolate import NearestNDInterpolator -from collections import defaultdict +from scipy.sparse import diags, linalg, lil_matrix +from lib.utils import setup_logger +log_file = snakemake.log[0] if snakemake.log else None +logger = setup_logger(log_file) -def find_boundary_vertices(vertices, faces): - """ - Find boundary vertices of a 3D mesh. - Args: - vertices (list of tuples): List of 3D points (x, y, z). - faces (list of tuples): List of triangular faces, where each face - is a tuple of three vertex indices (v1, v2, v3). - - Returns: - list: List of vertex indices that are boundary vertices, sorted in ascending order. +def get_terminal_indices(sdt, min_dist, max_dist, min_vertices, boundary_mask): """ - edge_count = defaultdict(int) - # Step 1: Count edge occurrences - for face in faces: - # Extract edges from the face, ensure consistent ordering (min, max) - edges = [ - tuple(sorted((face[0], face[1]))), - tuple(sorted((face[1], face[2]))), - tuple(sorted((face[2], face[0]))), - ] - for edge in edges: - edge_count[edge] += 1 - # Step 2: Identify boundary edges - boundary_edges = [edge for edge, count in edge_count.items() if count == 1] - # Step 3: Collect boundary vertices - boundary_vertices = set() - for edge in boundary_edges: - boundary_vertices.update(edge) - # Convert the set to a sorted list (array) - return sorted(boundary_vertices) + 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 + swept from min_dist to max_dist, until the min_vertices is achieved, else an + exception is thrown.""" + + for dist in np.linspace(min_dist, max_dist, 20): + indices = np.where((sdt < dist) & (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 {max_dist}mm of the terminal mask" + ) def solve_laplace_beltrami_open_mesh(vertices, faces, boundary_conditions=None): @@ -50,8 +35,11 @@ def solve_laplace_beltrami_open_mesh(vertices, faces, boundary_conditions=None): Returns: solution (np.ndarray): Array of shape (n_vertices,) with the solution values. """ + logger.info("solve_laplace_beltrami_open_mesh") n_vertices = vertices.shape[0] + logger.info(f"n_vertices: {n_vertices}") # Step 1: Compute cotangent weights + logger.info("Computing cotangent weights") weights = lil_matrix((n_vertices, n_vertices)) for tri in faces: v0, v1, v2 = vertices[tri[0]], vertices[tri[1]], vertices[tri[2]] @@ -75,8 +63,10 @@ def solve_laplace_beltrami_open_mesh(vertices, faces, boundary_conditions=None): weights[tri[2], tri[1]] += cot0 / 2 weights[tri[2], tri[0]] += cot1 / 2 weights[tri[0], tri[2]] += cot1 / 2 + logger.info("weights.tocsr()") weights = weights.tocsr() # Step 2: Handle boundaries for open meshes + logger.info("Handle boundaries for open meshes") diagonal = weights.sum(axis=1).A1 # Ensure no zero entries in diagonal to avoid singular matrix issues diagonal[diagonal < 1e-12] = 1e-12 @@ -86,12 +76,14 @@ def solve_laplace_beltrami_open_mesh(vertices, faces, boundary_conditions=None): boundary_indices = list(boundary_conditions.keys()) boundary_values = np.array(list(boundary_conditions.values())) free_indices = np.setdiff1d(np.arange(n_vertices), boundary_indices) + logger.info("Setting boundary conditions") b = np.zeros(n_vertices) for idx, value in boundary_conditions.items(): laplacian[idx, :] = 0 laplacian[idx, idx] = 1 b[idx] = value # Step 3: Solve the Laplace-Beltrami equation + logger.info("Solve the Laplace-Beltrami equation") solution = np.zeros(n_vertices) if len(free_indices) > 0: free_laplacian = laplacian[free_indices][:, free_indices] @@ -101,67 +93,57 @@ def solve_laplace_beltrami_open_mesh(vertices, faces, boundary_conditions=None): ) solution[boundary_indices] = boundary_values try: + logger.info("about to solve") solution[free_indices] = linalg.spsolve(free_laplacian, free_b) + logger.info("done solve") except linalg.MatrixRankWarning: - print("Warning: Laplacian matrix is singular or ill-conditioned.") + logger.info("Warning: Laplacian matrix is singular or ill-conditioned.") solution[free_indices] = np.zeros(len(free_indices)) else: solution[boundary_indices] = boundary_values return solution +logger.info( + "Loading in surface, boundary mask, and src/sink signed distance transforms" +) + surf = nib.load(snakemake.input.surf_gii) vertices = surf.agg_data("NIFTI_INTENT_POINTSET") faces = surf.agg_data("NIFTI_INTENT_TRIANGLE") +boundary_mask = nib.load(snakemake.input.boundary).agg_data() +src_sdt = nib.load(snakemake.input.src_sdt).agg_data() +sink_sdt = nib.load(snakemake.input.sink_sdt).agg_data() -# get source/sink vertices by nearest neighbour (only of edge vertices) -boundary_vertices = np.array(find_boundary_vertices(vertices, faces)) -seg = nib.load(snakemake.input.seg) - -src_AP = np.array( - np.where(seg.get_fdata() == snakemake.params.srcsink_labels["AP"]["src"]) -) -sink_AP = np.array( - np.where(seg.get_fdata() == snakemake.params.srcsink_labels["AP"]["sink"]) -) -src_PD = np.array( - np.where(seg.get_fdata() == snakemake.params.srcsink_labels["PD"]["src"]) -) -sink_PD = np.array( - np.where(seg.get_fdata() == snakemake.params.srcsink_labels["PD"]["sink"]) +src_indices = get_terminal_indices( + src_sdt, + snakemake.params.min_dist_threshold, + snakemake.params.max_dist_threshold, + snakemake.params.min_terminal_vertices, + boundary_mask, ) -# apply affine -src_AP = (seg.affine @ np.vstack([src_AP, np.ones([1, src_AP.shape[1]])]))[:3, :] -sink_AP = (seg.affine @ np.vstack([sink_AP, np.ones([1, sink_AP.shape[1]])]))[:3, :] -src_PD = (seg.affine @ np.vstack([src_PD, np.ones([1, src_PD.shape[1]])]))[:3, :] -sink_PD = (seg.affine @ np.vstack([sink_PD, np.ones([1, sink_PD.shape[1]])]))[:3, :] - -vals = np.hstack( - [ - np.ones([src_AP.shape[1]]) * 10, - np.ones([sink_AP.shape[1]]) * 11, - np.ones([src_PD.shape[1]]) * 20, - np.ones([sink_PD.shape[1]]) * 21, - ] +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, ) -interpol = NearestNDInterpolator(np.hstack([src_AP, sink_AP, src_PD, sink_PD]).T, vals) -boundary_values = np.array(interpol(vertices[boundary_vertices, :])).astype(int) -APinds = np.array(np.where(boundary_values < 12)[0]).astype(int) -boundary_conditions = dict(zip(boundary_vertices[APinds], boundary_values[APinds] - 10)) -APcoords = solve_laplace_beltrami_open_mesh(vertices, faces, boundary_conditions) -PDinds = np.array(np.where(boundary_values > 12)[0]).astype(int) -boundary_conditions = dict(zip(boundary_vertices[PDinds], boundary_values[PDinds] - 20)) -PDcoords = solve_laplace_beltrami_open_mesh(vertices, faces, boundary_conditions) +logger.info(f"# of sink boundary vertices: {len(sink_indices)}") -data_array = nib.gifti.GiftiDataArray(data=APcoords.astype(np.float32)) -image = nib.gifti.GiftiImage() -image.add_gifti_data_array(data_array) -nib.save(image, snakemake.output.coords_AP) +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)) + + +coords = solve_laplace_beltrami_open_mesh(vertices, faces, boundary_conditions) -data_array = nib.gifti.GiftiDataArray(data=PDcoords.astype(np.float32)) +data_array = nib.gifti.GiftiDataArray(data=coords.astype(np.float32)) image = nib.gifti.GiftiImage() image.add_gifti_data_array(data_array) -nib.save(image, snakemake.output.coords_PD) +nib.save(image, snakemake.output.coords) From 492817cdc75afee0a299d471cae409fa6b8505b2 Mon Sep 17 00:00:00 2001 From: Jordan DeKraker Date: Mon, 10 Feb 2025 16:54:43 -0500 Subject: [PATCH 2/5] merged dev into template-v2.0.0 --- hippunfold/workflow/rules/coords.smk | 45 +++--- hippunfold/workflow/rules/native_surf.smk | 119 ++++++++++++++-- .../workflow/scripts/get_boundary_vertices.py | 97 +++++++++++++ .../workflow/scripts/laplace_beltrami.py | 132 ++++++++---------- 4 files changed, 285 insertions(+), 108 deletions(-) create mode 100644 hippunfold/workflow/scripts/get_boundary_vertices.py diff --git a/hippunfold/workflow/rules/coords.smk b/hippunfold/workflow/rules/coords.smk index d585166a..e93711d9 100644 --- a/hippunfold/workflow/rules/coords.smk +++ b/hippunfold/workflow/rules/coords.smk @@ -24,21 +24,13 @@ def get_gm_labels(wildcards): return lbl_list -def get_sink_labels(wildcards): +def get_src_sink_labels(wildcards): lbl_list = " ".join( [ str(lbl) - for lbl in config["laplace_labels"][wildcards.label][wildcards.dir]["sink"] - ] - ) - return lbl_list - - -def get_src_labels(wildcards): - lbl_list = " ".join( - [ - str(lbl) - for lbl in config["laplace_labels"][wildcards.label][wildcards.dir]["src"] + for lbl in config["laplace_labels"][wildcards.label][wildcards.dir][ + wildcards.srcsink + ] ] ) return lbl_list @@ -103,11 +95,11 @@ def get_inputs_laplace(wildcards): return files -rule get_sink_mask: +rule get_src_sink_mask: input: labelmap=get_labels_for_laplace, params: - labels=get_sink_labels, + labels=get_src_sink_labels, output: mask=bids( root=work, @@ -115,7 +107,7 @@ rule get_sink_mask: suffix="mask.nii.gz", space="corobl", dir="{dir}", - desc="sink", + desc="{srcsink,src|sink}", hemi="{hemi}", label="{label}", **inputs.subj_wildcards, @@ -130,19 +122,28 @@ rule get_sink_mask: "c3d {input} -background -1 -retain-labels {params} -binarize {output}" -rule get_src_mask: +rule get_src_sink_sdt: + """calculate signed distance transform (negative inside, positive outside)""" input: - labelmap=get_labels_for_laplace, - params: - labels=get_src_labels, - output: mask=bids( root=work, datatype="coords", suffix="mask.nii.gz", space="corobl", dir="{dir}", - desc="src", + desc="{srcsink}", + hemi="{hemi}", + label="{label}", + **inputs.subj_wildcards, + ), + output: + sdt=bids( + root=work, + datatype="coords", + suffix="sdt.nii.gz", + space="corobl", + dir="{dir}", + desc="{srcsink}", hemi="{hemi}", label="{label}", **inputs.subj_wildcards, @@ -154,7 +155,7 @@ rule get_src_mask: group: "subj" shell: - "c3d {input} -background -1 -retain-labels {params} -binarize {output}" + "c3d {input} -sdt -o {output}" rule get_nan_mask: diff --git a/hippunfold/workflow/rules/native_surf.smk b/hippunfold/workflow/rules/native_surf.smk index 793d7dbb..3e691a9d 100644 --- a/hippunfold/workflow/rules/native_surf.smk +++ b/hippunfold/workflow/rules/native_surf.smk @@ -66,7 +66,7 @@ rule gen_native_mesh: params: threshold=lambda wildcards: surf_thresholds[wildcards.surfname], decimate_opts={ - "reduction": 0.5, + "reduction": 0.9, "feature_angle": 25, "preserve_topology": True, }, @@ -173,7 +173,7 @@ rule smooth_surface: # --- creating unfold surface from native anatomical, including post-processing -rule laplace_beltrami: +rule get_boundary_vertices: input: surf_gii=bids( root=root, @@ -184,25 +184,122 @@ rule laplace_beltrami: label="{label}", **inputs.subj_wildcards, ), - seg=get_labels_for_laplace, - params: - srcsink_labels=lambda wildcards: config["laplace_labels"][wildcards.label], output: - coords_AP=bids( + label_gii=bids( + root=work, + datatype="surf", + suffix="boundary.label.gii", + space="corobl", + hemi="{hemi}", + label="{label}", + **inputs.subj_wildcards, + ), + group: + "subj" + container: + config["singularity"]["autotop"] + conda: + "../envs/pyvista.yaml" + script: + "../scripts/get_boundary_vertices.py" + + +rule map_src_sink_sdt_to_surf: + """ Maps the distance to src/sink mask """ + input: + surf_gii=bids( + root=root, + datatype="surf", + suffix="midthickness.surf.gii", + space="corobl", + hemi="{hemi}", + label="{label}", + **inputs.subj_wildcards, + ), + sdt=bids( root=work, datatype="coords", - dir="AP", - suffix="coords.shape.gii", - desc="laplace", + suffix="sdt.nii.gz", space="corobl", + dir="{dir}", + desc="{srcsink}", hemi="{hemi}", label="{label}", **inputs.subj_wildcards, ), - coords_PD=bids( + output: + sdt=bids( + root=work, + datatype="surf", + suffix="sdt.shape.gii", + space="corobl", + hemi="{hemi}", + dir="{dir}", + desc="{srcsink}", + label="{label}", + **inputs.subj_wildcards, + ), + container: + config["singularity"]["autotop"] + conda: + "../envs/workbench.yaml" + group: + "subj" + shell: + "wb_command -volume-to-surface-mapping {input.sdt} {input.surf_gii} {output.sdt} -trilinear" + + +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_sdt=bids( + root=work, + datatype="surf", + suffix="sdt.shape.gii", + space="corobl", + hemi="{hemi}", + dir="{dir}", + desc="src", + label="{label}", + **inputs.subj_wildcards, + ), + sink_sdt=bids( + root=work, + datatype="surf", + suffix="sdt.shape.gii", + space="corobl", + hemi="{hemi}", + dir="{dir}", + desc="sink", + label="{label}", + **inputs.subj_wildcards, + ), + boundary=bids( + root=work, + datatype="surf", + suffix="boundary.label.gii", + space="corobl", + hemi="{hemi}", + label="{label}", + **inputs.subj_wildcards, + ), + params: + min_dist_threshold=0.3, + max_dist_threshold=1, + min_terminal_vertices=5, + output: + coords=bids( root=work, datatype="coords", - dir="PD", + dir="{dir}", suffix="coords.shape.gii", desc="laplace", space="corobl", diff --git a/hippunfold/workflow/scripts/get_boundary_vertices.py b/hippunfold/workflow/scripts/get_boundary_vertices.py new file mode 100644 index 00000000..1e5101d3 --- /dev/null +++ b/hippunfold/workflow/scripts/get_boundary_vertices.py @@ -0,0 +1,97 @@ +import pyvista as pv +import numpy as np +import nibabel as nib +import nibabel.gifti as gifti +from collections import defaultdict +from lib.utils import setup_logger + +# Setup logger +log_file = snakemake.log[0] if snakemake.log else None +logger = setup_logger(log_file) + + +def find_boundary_vertices(mesh): + """ + Find boundary vertices of a 3D mesh. + + Args: + mesh + Returns: + list: List of vertex indices that are boundary vertices, sorted in ascending order. + """ + vertices = mesh.points + faces = mesh.faces.reshape((-1, 4))[:, 1:4] # Extract triangle indices + + edge_count = defaultdict(int) + # Step 1: Count edge occurrences + for face in faces: + # Extract edges from the face, ensure consistent ordering (min, max) + edges = [ + tuple(sorted((face[0], face[1]))), + tuple(sorted((face[1], face[2]))), + tuple(sorted((face[2], face[0]))), + ] + for edge in edges: + edge_count[edge] += 1 + # Step 2: Identify boundary edges + boundary_edges = [edge for edge, count in edge_count.items() if count == 1] + # Step 3: Collect boundary vertices + boundary_vertices = set() + for edge in boundary_edges: + boundary_vertices.update(edge) + # Convert the set to a sorted list (array) + return np.array(sorted(boundary_vertices), dtype=np.int32) + + +def read_surface_from_gifti(surf_gii): + """Load a surface mesh from a GIFTI file.""" + surf = nib.load(surf_gii) + vertices = surf.agg_data("NIFTI_INTENT_POINTSET") + 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) + + +logger.info("Loading surface from GIFTI...") +surface = read_surface_from_gifti(snakemake.input.surf_gii) +logger.info(f"Surface loaded: {surface.n_points} vertices, {surface.n_faces} faces.") + + +logger.info("Find boundary vertices") +boundary_indices = find_boundary_vertices(surface) + +boundary_scalars = np.zeros(surface.n_points, dtype=np.int32) # Default is 0 +boundary_scalars[boundary_indices] = 1 # Set boundary vertices to 1 +logger.info( + f"Boundary scalar array created. {np.sum(boundary_scalars)} boundary vertices marked." +) + +logger.info("Saving GIFTI label file...") + +# Create a GIFTI label data array +gii_data = gifti.GiftiDataArray(boundary_scalars, intent="NIFTI_INTENT_LABEL") + +# Create a Label Table (LUT) +label_table = gifti.GiftiLabelTable() + +# Define Background label (key 0) +background_label = gifti.GiftiLabel( + key=0, red=1.0, green=1.0, blue=1.0, alpha=0.0 +) # Transparent +background_label.label = "Background" +label_table.labels.append(background_label) + +# Define Boundary label (key 1) +boundary_label = gifti.GiftiLabel( + key=1, red=1.0, green=0.0, blue=0.0, alpha=1.0 +) # Red color +boundary_label.label = "Boundary" +label_table.labels.append(boundary_label) + +# Assign label table to GIFTI image +gii_img = gifti.GiftiImage(darrays=[gii_data], labeltable=label_table) + +# Save the label file +gii_img.to_filename(snakemake.output.label_gii) +logger.info(f"GIFTI label file saved as '{snakemake.output.label_gii}'.") diff --git a/hippunfold/workflow/scripts/laplace_beltrami.py b/hippunfold/workflow/scripts/laplace_beltrami.py index de687ff1..e9cfabc8 100644 --- a/hippunfold/workflow/scripts/laplace_beltrami.py +++ b/hippunfold/workflow/scripts/laplace_beltrami.py @@ -1,41 +1,26 @@ import numpy as np -from scipy.sparse import diags, linalg, lil_matrix import nibabel as nib -from scipy.interpolate import NearestNDInterpolator -from collections import defaultdict +from scipy.sparse import diags, linalg, lil_matrix +from lib.utils import setup_logger +log_file = snakemake.log[0] if snakemake.log else None +logger = setup_logger(log_file) -def find_boundary_vertices(vertices, faces): - """ - Find boundary vertices of a 3D mesh. - Args: - vertices (list of tuples): List of 3D points (x, y, z). - faces (list of tuples): List of triangular faces, where each face - is a tuple of three vertex indices (v1, v2, v3). - - Returns: - list: List of vertex indices that are boundary vertices, sorted in ascending order. +def get_terminal_indices(sdt, min_dist, max_dist, min_vertices, boundary_mask): """ - edge_count = defaultdict(int) - # Step 1: Count edge occurrences - for face in faces: - # Extract edges from the face, ensure consistent ordering (min, max) - edges = [ - tuple(sorted((face[0], face[1]))), - tuple(sorted((face[1], face[2]))), - tuple(sorted((face[2], face[0]))), - ] - for edge in edges: - edge_count[edge] += 1 - # Step 2: Identify boundary edges - boundary_edges = [edge for edge, count in edge_count.items() if count == 1] - # Step 3: Collect boundary vertices - boundary_vertices = set() - for edge in boundary_edges: - boundary_vertices.update(edge) - # Convert the set to a sorted list (array) - return sorted(boundary_vertices) + 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 + swept from min_dist to max_dist, until the min_vertices is achieved, else an + exception is thrown.""" + + for dist in np.linspace(min_dist, max_dist, 20): + indices = np.where((sdt < dist) & (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 {max_dist}mm of the terminal mask" + ) def solve_laplace_beltrami_open_mesh(vertices, faces, boundary_conditions=None): @@ -50,8 +35,11 @@ def solve_laplace_beltrami_open_mesh(vertices, faces, boundary_conditions=None): Returns: solution (np.ndarray): Array of shape (n_vertices,) with the solution values. """ + logger.info("solve_laplace_beltrami_open_mesh") n_vertices = vertices.shape[0] + logger.info(f"n_vertices: {n_vertices}") # Step 1: Compute cotangent weights + logger.info("Computing cotangent weights") weights = lil_matrix((n_vertices, n_vertices)) for tri in faces: v0, v1, v2 = vertices[tri[0]], vertices[tri[1]], vertices[tri[2]] @@ -75,8 +63,10 @@ def solve_laplace_beltrami_open_mesh(vertices, faces, boundary_conditions=None): weights[tri[2], tri[1]] += cot0 / 2 weights[tri[2], tri[0]] += cot1 / 2 weights[tri[0], tri[2]] += cot1 / 2 + logger.info("weights.tocsr()") weights = weights.tocsr() # Step 2: Handle boundaries for open meshes + logger.info("Handle boundaries for open meshes") diagonal = weights.sum(axis=1).A1 # Ensure no zero entries in diagonal to avoid singular matrix issues diagonal[diagonal < 1e-12] = 1e-12 @@ -86,12 +76,14 @@ def solve_laplace_beltrami_open_mesh(vertices, faces, boundary_conditions=None): boundary_indices = list(boundary_conditions.keys()) boundary_values = np.array(list(boundary_conditions.values())) free_indices = np.setdiff1d(np.arange(n_vertices), boundary_indices) + logger.info("Setting boundary conditions") b = np.zeros(n_vertices) for idx, value in boundary_conditions.items(): laplacian[idx, :] = 0 laplacian[idx, idx] = 1 b[idx] = value # Step 3: Solve the Laplace-Beltrami equation + logger.info("Solve the Laplace-Beltrami equation") solution = np.zeros(n_vertices) if len(free_indices) > 0: free_laplacian = laplacian[free_indices][:, free_indices] @@ -101,67 +93,57 @@ def solve_laplace_beltrami_open_mesh(vertices, faces, boundary_conditions=None): ) solution[boundary_indices] = boundary_values try: + logger.info("about to solve") solution[free_indices] = linalg.spsolve(free_laplacian, free_b) + logger.info("done solve") except linalg.MatrixRankWarning: - print("Warning: Laplacian matrix is singular or ill-conditioned.") + logger.info("Warning: Laplacian matrix is singular or ill-conditioned.") solution[free_indices] = np.zeros(len(free_indices)) else: solution[boundary_indices] = boundary_values return solution +logger.info( + "Loading in surface, boundary mask, and src/sink signed distance transforms" +) + surf = nib.load(snakemake.input.surf_gii) vertices = surf.agg_data("NIFTI_INTENT_POINTSET") faces = surf.agg_data("NIFTI_INTENT_TRIANGLE") +boundary_mask = nib.load(snakemake.input.boundary).agg_data() +src_sdt = nib.load(snakemake.input.src_sdt).agg_data() +sink_sdt = nib.load(snakemake.input.sink_sdt).agg_data() -# get source/sink vertices by nearest neighbour (only of edge vertices) -boundary_vertices = np.array(find_boundary_vertices(vertices, faces)) -seg = nib.load(snakemake.input.seg) - -src_AP = np.array( - np.where(seg.get_fdata() == snakemake.params.srcsink_labels["AP"]["src"]) -) -sink_AP = np.array( - np.where(seg.get_fdata() == snakemake.params.srcsink_labels["AP"]["sink"]) -) -src_PD = np.array( - np.where(seg.get_fdata() == snakemake.params.srcsink_labels["PD"]["src"]) -) -sink_PD = np.array( - np.where(seg.get_fdata() == snakemake.params.srcsink_labels["PD"]["sink"]) +src_indices = get_terminal_indices( + src_sdt, + snakemake.params.min_dist_threshold, + snakemake.params.max_dist_threshold, + snakemake.params.min_terminal_vertices, + boundary_mask, ) -# apply affine -src_AP = (seg.affine @ np.vstack([src_AP, np.ones([1, src_AP.shape[1]])]))[:3, :] -sink_AP = (seg.affine @ np.vstack([sink_AP, np.ones([1, sink_AP.shape[1]])]))[:3, :] -src_PD = (seg.affine @ np.vstack([src_PD, np.ones([1, src_PD.shape[1]])]))[:3, :] -sink_PD = (seg.affine @ np.vstack([sink_PD, np.ones([1, sink_PD.shape[1]])]))[:3, :] - -vals = np.hstack( - [ - np.ones([src_AP.shape[1]]) * 10, - np.ones([sink_AP.shape[1]]) * 11, - np.ones([src_PD.shape[1]]) * 20, - np.ones([sink_PD.shape[1]]) * 21, - ] +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, ) -interpol = NearestNDInterpolator(np.hstack([src_AP, sink_AP, src_PD, sink_PD]).T, vals) -boundary_values = np.array(interpol(vertices[boundary_vertices, :])).astype(int) -APinds = np.array(np.where(boundary_values < 12)[0]).astype(int) -boundary_conditions = dict(zip(boundary_vertices[APinds], boundary_values[APinds] - 10)) -APcoords = solve_laplace_beltrami_open_mesh(vertices, faces, boundary_conditions) -PDinds = np.array(np.where(boundary_values > 12)[0]).astype(int) -boundary_conditions = dict(zip(boundary_vertices[PDinds], boundary_values[PDinds] - 20)) -PDcoords = solve_laplace_beltrami_open_mesh(vertices, faces, boundary_conditions) +logger.info(f"# of sink boundary vertices: {len(sink_indices)}") -data_array = nib.gifti.GiftiDataArray(data=APcoords.astype(np.float32)) -image = nib.gifti.GiftiImage() -image.add_gifti_data_array(data_array) -nib.save(image, snakemake.output.coords_AP) +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)) + + +coords = solve_laplace_beltrami_open_mesh(vertices, faces, boundary_conditions) -data_array = nib.gifti.GiftiDataArray(data=PDcoords.astype(np.float32)) +data_array = nib.gifti.GiftiDataArray(data=coords.astype(np.float32)) image = nib.gifti.GiftiImage() image.add_gifti_data_array(data_array) -nib.save(image, snakemake.output.coords_PD) +nib.save(image, snakemake.output.coords) From 16b31783722a2145cbdc498959894ec0249af153 Mon Sep 17 00:00:00 2001 From: Jordan DeKraker Date: Mon, 10 Feb 2025 17:01:52 -0500 Subject: [PATCH 3/5] Update python-testing.yml for dsegtissue --- .github/workflows/python-testing.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/python-testing.yml b/.github/workflows/python-testing.yml index 68a3c26e..40da2eb1 100644 --- a/.github/workflows/python-testing.yml +++ b/.github/workflows/python-testing.yml @@ -96,11 +96,11 @@ jobs: - name: Test manualseg bids, with path override run: | - poetry run hippunfold test_data/bids_manualseg test_out participant -np --modality manualseg --path_manualseg test_data/bids_manualseg/sub-{subject}_hemi-{hemi}_dseg.nii.gz + poetry run hippunfold test_data/bids_manualseg test_out participant -np --modality dsegtissue --path_dsegtissue test_data/bids_manualseg/sub-{subject}_hemi-{hemi}_dseg.nii.gz - name: Test manualseg bids, one hemisphere hemi run: | - poetry run hippunfold test_data/bids_manualseg_1hemi test_out participant -np --modality manualseg --hemi L + poetry run hippunfold test_data/bids_manualseg_1hemi test_out participant -np --modality dsegtissue --hemi L - name: Test T2w with T1w template registration run: | From 0cee11a1bbfdb67cae3a8ffc9012f65fb5353c4a Mon Sep 17 00:00:00 2001 From: Jordan DeKraker Date: Mon, 10 Feb 2025 17:03:32 -0500 Subject: [PATCH 4/5] dsegtissue in dryruns to current branch --- .github/workflows/python-testing.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/python-testing.yml b/.github/workflows/python-testing.yml index 68a3c26e..40da2eb1 100644 --- a/.github/workflows/python-testing.yml +++ b/.github/workflows/python-testing.yml @@ -96,11 +96,11 @@ jobs: - name: Test manualseg bids, with path override run: | - poetry run hippunfold test_data/bids_manualseg test_out participant -np --modality manualseg --path_manualseg test_data/bids_manualseg/sub-{subject}_hemi-{hemi}_dseg.nii.gz + poetry run hippunfold test_data/bids_manualseg test_out participant -np --modality dsegtissue --path_dsegtissue test_data/bids_manualseg/sub-{subject}_hemi-{hemi}_dseg.nii.gz - name: Test manualseg bids, one hemisphere hemi run: | - poetry run hippunfold test_data/bids_manualseg_1hemi test_out participant -np --modality manualseg --hemi L + poetry run hippunfold test_data/bids_manualseg_1hemi test_out participant -np --modality dsegtissue --hemi L - name: Test T2w with T1w template registration run: | From e4c1562294f8ae8bb605cfdc1ff24068623ae095 Mon Sep 17 00:00:00 2001 From: Ali Khan Date: Mon, 10 Feb 2025 20:02:31 -0500 Subject: [PATCH 5/5] another rename manualseg in tests --- .github/workflows/python-testing.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/python-testing.yml b/.github/workflows/python-testing.yml index 40da2eb1..7b4389a3 100644 --- a/.github/workflows/python-testing.yml +++ b/.github/workflows/python-testing.yml @@ -96,11 +96,11 @@ jobs: - name: Test manualseg bids, with path override run: | - poetry run hippunfold test_data/bids_manualseg test_out participant -np --modality dsegtissue --path_dsegtissue test_data/bids_manualseg/sub-{subject}_hemi-{hemi}_dseg.nii.gz + poetry run hippunfold test_data/bids_dsegtissue test_out participant -np --modality dsegtissue --path_dsegtissue test_data/bids_dsegtissue/sub-{subject}_hemi-{hemi}_dseg.nii.gz - name: Test manualseg bids, one hemisphere hemi run: | - poetry run hippunfold test_data/bids_manualseg_1hemi test_out participant -np --modality dsegtissue --hemi L + poetry run hippunfold test_data/bids_dsegtissue_1hemi test_out participant -np --modality dsegtissue --hemi L - name: Test T2w with T1w template registration run: |