Skip to content

Commit

Permalink
lint
Browse files Browse the repository at this point in the history
  • Loading branch information
Jordan DeKraker committed Jan 21, 2025
1 parent 5d98337 commit fff3695
Show file tree
Hide file tree
Showing 2 changed files with 31 additions and 17 deletions.
38 changes: 25 additions & 13 deletions hippunfold/workflow/scripts/laplace-beltrami.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,38 +5,46 @@
from collections import defaultdict
import networkx as nx


def largest_connected_component(vertices: np.ndarray, faces: np.ndarray):
"""
Returns the vertices, faces, and original indices of the largest connected component of a 3D surface mesh.
Parameters:
vertices (np.ndarray): Array of vertex coordinates (N x 3).
faces (np.ndarray): Array of face indices (M x 3).
Returns:
tuple: (new_vertices, new_faces, largest_component), where new_vertices are the vertex coordinates of the largest component,
new_faces are the face indices adjusted to the new vertex order, and largest_component contains the indices of the original vertices.
"""
# Build adjacency graph from face connectivity
G = nx.Graph()
for face in faces:
G.add_edges_from([(face[i], face[j]) for i in range(3) for j in range(i + 1, 3)])

G.add_edges_from(
[(face[i], face[j]) for i in range(3) for j in range(i + 1, 3)]
)

# Find connected components
components = list(nx.connected_components(G))

# Select the largest connected component
largest_component = max(components, key=len)
largest_component = np.array(list(largest_component))

# Create a mapping from old vertex indices to new ones
index_map = {old_idx: new_idx for new_idx, old_idx in enumerate(largest_component)}

# Filter the vertices and faces
new_vertices = vertices[largest_component]
new_faces = np.array([[index_map[v] for v in face] for face in faces
if all(v in index_map for v in face)])

new_faces = np.array(
[
[index_map[v] for v in face]
for face in faces
if all(v in index_map for v in face)
]
)

return new_vertices, new_faces, largest_component


Expand Down Expand Up @@ -149,7 +157,7 @@ def solve_laplace_beltrami_open_mesh(vertices, faces, boundary_conditions=None):
vertices_orig = surf.agg_data("NIFTI_INTENT_POINTSET")
faces = surf.agg_data("NIFTI_INTENT_TRIANGLE")

vertices, faces, i_concomp = largest_connected_component(vertices_orig,faces)
vertices, faces, i_concomp = largest_connected_component(vertices_orig, faces)


# get source/sink vertices by nearest neighbour (only of edge vertices)
Expand Down Expand Up @@ -185,11 +193,15 @@ def solve_laplace_beltrami_open_mesh(vertices, faces, boundary_conditions=None):
APinds = np.array(np.where(boundary_values < 12)[0]).astype(int)
boundary_conditions = dict(zip(boundary_vertices[APinds], boundary_values[APinds] - 10))
APcoords = np.nans((len(vertices_orig)))
APcoords[i_concomp] = solve_laplace_beltrami_open_mesh(vertices, faces, boundary_conditions)
APcoords[i_concomp] = 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 = np.nans((len(vertices_orig)))
PDcoords[i_concomp] = solve_laplace_beltrami_open_mesh(vertices, faces, boundary_conditions)
PDcoords[i_concomp] = solve_laplace_beltrami_open_mesh(
vertices, faces, boundary_conditions
)

data_array = nib.gifti.GiftiDataArray(data=APcoords.astype(np.float32))
image = nib.gifti.GiftiImage()
Expand Down
10 changes: 6 additions & 4 deletions hippunfold/workflow/scripts/rewrite_vertices_to_flat.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,10 +6,12 @@
ap = nib.load(snakemake.input.coords_AP).darrays[0].data
pd = nib.load(snakemake.input.coords_PD).darrays[0].data

vertices[:, 0] = ap * -float(snakemake.params.vertspace["extent"][0]
) - float(snakemake.params.vertspace["origin"][0])
vertices[:, 1] = pd * float(snakemake.params.vertspace["extent"][1]
) - float(snakemake.params.vertspace["origin"][1])
vertices[:, 0] = ap * -float(snakemake.params.vertspace["extent"][0]) - float(
snakemake.params.vertspace["origin"][0]
)
vertices[:, 1] = pd * float(snakemake.params.vertspace["extent"][1]) - float(
snakemake.params.vertspace["origin"][1]
)
vertices[:, 2] = snakemake.params.z_level

nib.save(gii, snakemake.output.surf_gii)

0 comments on commit fff3695

Please sign in to comment.