From cae0f4991857702163d2e93c3382995c32aa285a Mon Sep 17 00:00:00 2001 From: Martin Reuter Date: Fri, 16 Jan 2026 18:15:59 +0100 Subject: [PATCH 1/8] add crtical points, boundary and improve level paths --- lapy/tria_mesh.py | 344 +++++++++++++++++++++++++++++++++++++--------- 1 file changed, 280 insertions(+), 64 deletions(-) diff --git a/lapy/tria_mesh.py b/lapy/tria_mesh.py index a40e5845..1d000546 100644 --- a/lapy/tria_mesh.py +++ b/lapy/tria_mesh.py @@ -329,6 +329,24 @@ def is_manifold(self): """ return np.max(self.adj_sym.data) <= 2 + def is_boundary(self) -> np.ndarray: + """Check which vertices are on the boundary. + + Returns + ------- + np.ndarray + Boolean array of shape (n_vertices,) where True indicates the vertex + is on a boundary edge (an edge that belongs to only one triangle). + """ + # Boundary edges have value 1 in the symmetric adjacency matrix + boundary_edges = self.adj_sym == 1 + # Get vertices that are part of any boundary edge + boundary_vertices = np.zeros(self.v.shape[0], dtype=bool) + rows, cols = boundary_edges.nonzero() + boundary_vertices[rows] = True + boundary_vertices[cols] = True + return boundary_vertices + def is_oriented(self) -> bool: """Check if triangle mesh is oriented. @@ -1014,6 +1032,112 @@ def curvature_tria( # find orthogonal umin next? return tumin2, tumax2, tcmin, tcmax + def critical_points(self, vfunc: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + """Compute critical points (extrema and saddles) of a function on mesh vertices. + + A minimum is a vertex where all neighbor values are larger. + A maximum is a vertex where all neighbor values are smaller. + A saddle is a vertex with at least two regions of neighbors with larger values, + and two with smaller values. Boundary vertices assume a virtual edge outside + the mesh that closes the boundary loop around the vertex. + + Parameters + ---------- + vfunc : np.ndarray + Real-valued function defined on mesh vertices, shape (n_vertices,). + + Returns + ------- + minima : np.ndarray + Array of vertex indices that are local minima, shape (n_minima,). + maxima : np.ndarray + Array of vertex indices that are local maxima, shape (n_maxima,). + saddles : np.ndarray + Array of vertex indices that are saddles (all orders), shape (n_saddles,). + saddle_orders : np.ndarray + Array of saddle orders for each saddle vertex (same length as saddles), shape (n_saddles,). + Order 2 = simple saddle (4 sign flips), order 3+ = higher-order saddles. + + Notes + ----- + The algorithm works by: + 1. For each vertex, compute difference: neighbor_value - vertex_value + 2. Minima: all differences are positive (all neighbors higher) + 3. Maxima: all differences are negative (all neighbors lower) + 4. Saddles: count sign flips across opposite edges in triangles at vertex + - Regular point: 2 sign flips + - Simple saddle (order 2): 4 sign flips + - Higher-order saddle (order n): 2n sign flips, order = n_flips / 2 + 5. Tie-breaker: when two vertices have equal function values, the vertex + with the higher vertex ID is treated as slightly larger to remove ambiguity. + + The implementation is fully vectorized using sparse matrices and numpy operations. + """ + vfunc = np.asarray(vfunc) + if len(vfunc) != self.v.shape[0]: + raise ValueError("vfunc length must match number of vertices") + n_vertices = self.v.shape[0] + + # Use SYMMETRIC adjacency matrix to get ALL edges (including boundary edges in both directions) + # COMPUTE EDGE SIGNS ONCE for all neighbor relationships + rows, cols = self.adj_sym.nonzero() + edge_diffs = vfunc[cols] - vfunc[rows] + edge_signs = np.sign(edge_diffs) + # Tie-breaker: when function values are equal, vertex with higher ID is larger + zero_mask = edge_signs == 0 + edge_signs[zero_mask] = np.sign(cols[zero_mask] - rows[zero_mask]) + # Create sparse matrix of edge signs for O(1) lookup + edge_sign_matrix = sparse.csr_matrix( + (edge_signs, (rows, cols)), + shape=(n_vertices, n_vertices) + ) + + # CLASSIFY MINIMA AND MAXIMA + # Compute min and max edge sign per vertex (row-wise) + # Note: edge_sign_matrix only contains +1/-1 (never 0 due to tie-breaker) + # Implicit zeros in sparse matrix represent non-neighbors and can be ignored + min_signs = edge_sign_matrix.min(axis=1).toarray().flatten() + max_signs = edge_sign_matrix.max(axis=1).toarray().flatten() + # Minimum: all neighbor signs are positive (+1) + # If min in {0,1} and max=1, all actual neighbors are +1 (zeros are non-neighbors) + is_minimum = (min_signs > -1) & (max_signs == 1) + # Maximum: all neighbor signs are negative (-1) + # If min=-1 and max in {-1,0}, all actual neighbors are -1 (zeros are non-neighbors) + is_maximum = (min_signs == -1) & (max_signs < 1) + minima = np.where(is_minimum)[0] + maxima = np.where(is_maximum)[0] + + # COUNT SIGN FLIPS at opposite edge for saddle detection + sign_flips = np.zeros(n_vertices, dtype=int) + t0 = self.t[:, 0] + t1 = self.t[:, 1] + t2 = self.t[:, 2] + # For vertex 0, opposite edge is (v1, v2) + sign0_1 = np.array(edge_sign_matrix[t0, t1]).flatten() + sign0_2 = np.array(edge_sign_matrix[t0, t2]).flatten() + sign1_2 = np.array(edge_sign_matrix[t1, t2]).flatten() + flip0 = (sign0_1 * sign0_2) < 0 + np.add.at(sign_flips, t0[flip0], 1) + # For vertex 1, opposite edge is (v2, v0) + # sign(v1→v0) = -sign(v0→v1) = -sign0_1 + flip1 = (sign1_2 * (-sign0_1)) < 0 + np.add.at(sign_flips, t1[flip1], 1) + # For vertex 2, opposite edge is (v0, v1) + # sign(v2→v0) = -sign(v0→v2) = -sign0_2 + # sign(v2→v1) = -sign(v1→v2) = -sign1_2 + flip2 = (sign0_2 * sign1_2) < 0 + np.add.at(sign_flips, t2[flip2], 1) + + # CLASSIFY SADDLES + # Saddles have 4+ flips (regular points have 2, minima/maxima have 0) + # Boundary vertices can have 3 flips (or more) to be a saddle, assuming an additional flip outside + # Order = n_flips / 2 + is_saddle = sign_flips >= 3 + saddles = np.where(is_saddle)[0] + saddle_orders = (sign_flips[saddles] + 1) // 2 + + return minima, maxima, saddles, saddle_orders + def normalize_(self) -> None: """Normalize TriaMesh to unit surface area and centroid at the origin. @@ -1562,92 +1686,169 @@ def __reduce_edges_to_path( edges: np.ndarray, start_idx: Optional[int] = None, get_edge_idx: bool = False, - ) -> Union[np.ndarray, tuple[np.ndarray, np.ndarray]]: - """Reduce undirected unsorted edges (index pairs) to path (index array). + ) -> Union[list, tuple[list, list]]: + """Reduce undirected unsorted edges to ordered path(s). - Converts an unordered list of edge pairs into an ordered path by finding - a traversal from one endpoint to the other. + Converts unordered edge pairs into ordered paths by finding traversals. + Handles single open paths, closed loops, and multiple disconnected components. + Always returns lists to handle multiple components uniformly. Parameters ---------- edges : np.ndarray Array of shape (n, 2) with pairs of positive integer node indices - representing undirected edges. All indices from 0 to max(edges) must - be used and graph needs to be connected. Nodes cannot have more than - 2 neighbors. Requires exactly two endnodes (nodes with only one - neighbor). For closed loops, cut it open by removing one edge before - passing to this function. + representing undirected edges. Node indices can have gaps (i.e., not all + indices from 0 to max need to appear in the edges). Only nodes that + actually appear in edges are processed. start_idx : int or None, default=None - Node with only one neighbor to start path. If None, the node with - the smaller index will be selected as start point. + Preferred start node. If None, selects an endpoint (degree 1) automatically, + or arbitrary node for closed loops. Must be a node that appears in edges. get_edge_idx : bool, default=False - If True, also return index of edge into edges for each path segment. - The returned edge index list has length n, while path has length n+1. + If True, also return edge indices for each path segment. Returns ------- - path : np.ndarray - Array of length n+1 containing node indices as single path from start - to endpoint. - edge_idx : np.ndarray - Array of length n containing corresponding edge idx into edges for each - path segment. Only returned if get_edge_idx is True. + paths : list[np.ndarray] + List of ordered paths, one per connected component. For closed loops, + first node is repeated at end. + edge_idxs : list[np.ndarray] + List of edge index arrays, one per component. Only returned if + get_edge_idx=True. Raises ------ ValueError - If graph does not have exactly two endpoints. - If start_idx is not one of the endpoints. - If graph is not connected. + If start_idx is invalid. + If graph has degree > 2 (branching or self-intersections). """ - from scipy.sparse.csgraph import shortest_path - - # Extract node indices and create a sparse adjacency matrix edges = np.array(edges) + if edges.shape[0] == 0: + return ([], []) if get_edge_idx else [] + + # Build adjacency matrix i = np.column_stack((edges[:, 0], edges[:, 1])).reshape(-1) j = np.column_stack((edges[:, 1], edges[:, 0])).reshape(-1) dat = np.ones(i.shape) n = edges.max() + 1 adj_matrix = sparse.csr_matrix((dat, (i, j)), shape=(n, n)) - # Find the degree of each node, sum over rows to get degree degrees = np.asarray(adj_matrix.sum(axis=1)).ravel() - endpoints = np.where(degrees == 1)[0] - if len(endpoints) != 2: - raise ValueError( - "The graph does not have exactly two endpoints; invalid input." - ) - if not start_idx: - start_idx = endpoints[0] - else: - if not np.isin(start_idx, endpoints): + + # Find connected components + n_comp, labels = sparse.csgraph.connected_components(adj_matrix, directed=False) + + # Build edge index lookup matrix + edge_dat = np.arange(edges.shape[0]) + 1 + edge_dat = np.column_stack((edge_dat, edge_dat)).reshape(-1) + eidx_matrix = sparse.csr_matrix((edge_dat, (i, j)), shape=(n, n)) + + paths = [] + edge_idxs = [] + + for comp_id in range(n_comp): + comp_nodes = np.where(labels == comp_id)[0] + if len(comp_nodes) == 0: + continue + + comp_degrees = degrees[comp_nodes] + + # Skip isolated nodes (degree 0) that exist only due to matrix sizing. + # When edges use indices [0, 5, 10], we create a matrix of size 11x11. + # Indices [1,2,3,4,6,7,8,9] don't appear in any edges (have degree 0). + # connected_components treats each as a separate component, but they're + # not real nodes - just phantom entries from sizing matrix to max_index+1. + if np.all(comp_degrees == 0): + continue + + # SAFETY CHECK: Reject graphs with nodes of degree > 2 + # This single check catches all problematic cases: + # - Branching structures (Y-shape, star graph) + # - Self-intersections (figure-8, etc.) + # - Trees with multiple endpoints + # Valid graphs have only degree 1 (endpoints) and degree 2 (path nodes) + max_degree = comp_degrees.max() + if max_degree > 2: + high_degree_nodes = comp_nodes[comp_degrees > 2] raise ValueError( - f"start_idx {start_idx} must be one of the endpoints {endpoints}." + f"Component {comp_id}: found {len(high_degree_nodes)} node(s) with degree > 2 " + f"(max degree: {max_degree}). Degrees: {np.sort(comp_degrees)}. " + f"This indicates branching or self-intersecting structure. " + f"Only simple paths and simple closed loops are supported." ) - # Traverse the graph by computing shortest path - dist_matrix = shortest_path( - csgraph=adj_matrix, - directed=False, - indices=start_idx, - return_predecessors=False, - ) - if np.isinf(dist_matrix).any(): - raise ValueError( - "Ensure connected graph with indices from 0 to max_idx without gaps." - ) - # sort indices according to distance form start - path = dist_matrix.argsort() - # get edge idx of each segment from original list - enum = edges.shape[0] - dat = np.arange(enum) + 1 - dat = np.column_stack((dat, dat)).reshape(-1) - eidx_matrix = sparse.csr_matrix((dat, (i, j)), shape=(n, n)) - ei = path[0:-1] - ej = path[(np.arange(path.size - 1) + 1)] - eidx = np.asarray(eidx_matrix[ei, ej] - 1).ravel() + + # Determine if closed loop: all nodes have degree 2 + is_closed = np.all(comp_degrees == 2) + + # For closed loops: break one edge to convert to open path + if is_closed: + # Pick start node + if start_idx is not None and start_idx in comp_nodes: + start = start_idx + else: + start = comp_nodes[0] + + # Find neighbors of start node to break one edge + neighbors = adj_matrix[start, :].nonzero()[1] + neighbors_in_comp = [n for n in neighbors if n in comp_nodes] + + if len(neighbors_in_comp) < 2: + raise ValueError(f"Component {comp_id}: Closed loop node {start} should have 2 neighbors") + + # Break the edge from start to first neighbor (temporarily) + adj_matrix = adj_matrix.tolil() + break_neighbor = neighbors_in_comp[0] + adj_matrix[start, break_neighbor] = 0 + adj_matrix[break_neighbor, start] = 0 + adj_matrix = adj_matrix.tocsr() + + # Update degrees after breaking edge + degrees = np.asarray(adj_matrix.sum(axis=1)).ravel() + comp_degrees = degrees[comp_nodes] + + # Now handle as open path (both originally open and converted closed loops) + endpoints = comp_nodes[comp_degrees == 1] + + if len(endpoints) != 2: + raise ValueError( + f"Component {comp_id}: Expected 2 endpoints after breaking loop, found {len(endpoints)}" + ) + + # Select start node + if is_closed: + # For closed loops, start is already selected above + pass + elif start_idx is not None and start_idx in endpoints: + start = start_idx + else: + # For originally open paths, pick first endpoint + start = endpoints[0] + + # Use shortest_path to order nodes by distance from start + dist = sparse.csgraph.shortest_path(adj_matrix, indices=start, directed=False) + + if np.isinf(dist[comp_nodes]).any(): + raise ValueError(f"Component {comp_id} is not fully connected.") + + # Sort nodes by distance from start + path = comp_nodes[dist[comp_nodes].argsort()] + + # For closed loops: append start again to close the loop + if is_closed: + path = np.append(path, path[0]) + + paths.append(path) + + # Get edge indices if requested + if get_edge_idx: + ei = path[:-1] + ej = path[1:] + eidx = np.asarray(eidx_matrix[ei, ej] - 1).ravel() + edge_idxs.append(eidx) + + # Always return lists if get_edge_idx: - return path, eidx + return paths, edge_idxs else: - return path + return paths def level_path( @@ -1766,15 +1967,31 @@ def level_path( p2 = np.squeeze(p[edge_idxs[:, 1]]) llength = np.sqrt(((p1 - p2) ** 2).sum(1)).sum() # compute path from unordered, not-directed edge list - # and return path as list of points, and path length + # __reduce_edges_to_path now returns lists (can have multiple components) if get_tria_idx: - path, edge_idx = TriaMesh.__reduce_edges_to_path( - edge_idxs, get_edge_idx=get_tria_idx + paths, edge_idxs_list = TriaMesh.__reduce_edges_to_path( + edge_idxs, get_edge_idx=True ) + # level_path expects single component + if len(paths) != 1: + raise ValueError( + f"Found {len(paths)} disconnected level curves. " + f"Use extract_level_curves() for multiple components." + ) + path = paths[0] + edge_idx = edge_idxs_list[0] # translate local edge id to global tria id tria_idx = t_idx[edge_idx] else: - path = TriaMesh.__reduce_edges_to_path(edge_idxs, get_tria_idx) + paths = TriaMesh.__reduce_edges_to_path(edge_idxs, get_edge_idx=False) + # level_path expects single component + if len(paths) != 1: + raise ValueError( + f"Found {len(paths)} disconnected level curves. " + f"Use extract_level_curves() for multiple components." + ) + path = paths[0] + # remove duplicate vertices (happens when levelset hits a vertex almost # perfectly) path3d = p[path, :] @@ -1808,4 +2025,3 @@ def level_path( return path3d, llength, edges_vidxs, edges_relpos else: return path3d, llength - From d925bfc539e626a4451ec8693c0b2442f138658f Mon Sep 17 00:00:00 2001 From: Martin Reuter Date: Wed, 21 Jan 2026 17:49:17 +0100 Subject: [PATCH 2/8] allow to auto-determine if closed or not based on last/first point duplication when initialized --- lapy/polygon.py | 59 +++++++++++++++++++++++++++++++++++-------------- 1 file changed, 42 insertions(+), 17 deletions(-) diff --git a/lapy/polygon.py b/lapy/polygon.py index fbe10af9..2185a5d5 100644 --- a/lapy/polygon.py +++ b/lapy/polygon.py @@ -5,6 +5,7 @@ """ import logging import sys +from typing import Optional import numpy as np from scipy import sparse @@ -22,15 +23,21 @@ class Polygon: points : np.ndarray Array of shape (n, d) containing coordinates of polygon vertices in order, where d is 2 or 3 for 2D (x, y) or 3D (x, y, z) paths. - For closed polygons, the last point should not duplicate the first point. - closed : bool, default=False - If True, treats the path as a closed polygon. If False, treats it as - an open polyline. + If the last point duplicates the first point and closed is None, + the duplicate endpoint will be removed and the polygon will be + marked as closed automatically. + closed : bool or None, default=None + - If None (default): Auto-detect by checking if first and last points + are identical. If they are, removes duplicate endpoint and sets closed=True. + - If True: Treats the path as a closed polygon. If the last point duplicates + the first, it will be removed. + - If False: Treats the path as an open polyline regardless of endpoints. Attributes ---------- points : np.ndarray - Polygon vertex coordinates, shape (n_points, d). + Polygon vertex coordinates, shape (n_points, d). For closed polygons, + the last point does not duplicate the first. closed : bool Whether the polygon is closed or open. _is_2d : bool @@ -46,23 +53,22 @@ class Polygon: -------- >>> import numpy as np >>> from lapy.polygon import Polygon - >>> # Create a 2D closed polygon (square) - >>> square = np.array([[0, 0], [1, 0], [1, 1], [0, 1]]) - >>> poly = Polygon(square, closed=True) - >>> poly.is_2d() + >>> # Create a 2D closed polygon (square) - auto-detected + >>> square = np.array([[0, 0], [1, 0], [1, 1], [0, 1], [0, 0]]) + >>> poly = Polygon(square) # closed=None, will auto-detect and remove duplicate + >>> poly.closed True - >>> poly.length() - 4.0 - >>> # Create a 3D open path - >>> path_3d = np.array([[0, 0, 0], [1, 0, 0], [1, 1, 0], [1, 1, 1]]) - >>> poly3d = Polygon(path_3d, closed=False) - >>> poly3d.is_2d() + >>> poly.points.shape[0] + 4 + >>> # Explicitly open polygon + >>> path = np.array([[0, 0, 0], [1, 0, 0], [1, 1, 0], [1, 1, 1]]) + >>> poly3d = Polygon(path, closed=False) + >>> poly3d.closed False """ - def __init__(self, points: np.ndarray, closed: bool = False): + def __init__(self, points: np.ndarray, closed: Optional[bool] = None): self.points = np.array(points) - self.closed = closed # Validate non-empty polygon if self.points.size == 0: @@ -93,6 +99,25 @@ def __init__(self, points: np.ndarray, closed: bool = False): else: raise ValueError("Points should have 2 or 3 coordinates") + # Auto-detect closed polygons or handle explicit closed parameter + if closed is None: + # Auto-detect: check if first and last points are identical + if len(self.points) > 1 and np.allclose(self.points[0], self.points[-1]): + # Remove duplicate endpoint + self.points = self.points[:-1] + self.closed = True + else: + # Not closed (endpoints differ) + self.closed = False + elif closed: + # Explicitly closed: remove duplicate endpoint if present + if len(self.points) > 1 and np.allclose(self.points[0], self.points[-1]): + self.points = self.points[:-1] + self.closed = True + else: + # Explicitly open + self.closed = False + def is_2d(self) -> bool: """Check if the polygon is 2D. From a924c1b49ad3fb1b1c80d4a12106923280467d13 Mon Sep 17 00:00:00 2001 From: Martin Reuter Date: Wed, 21 Jan 2026 17:50:58 +0100 Subject: [PATCH 3/8] add extract_level_paths for multiple non-crossing/branching levelsets --- lapy/tria_mesh.py | 354 +++++++++++++++------- lapy/utils/tests/test_tria_mesh.py | 472 +++++++++++++++++++++++++++++ 2 files changed, 718 insertions(+), 108 deletions(-) diff --git a/lapy/tria_mesh.py b/lapy/tria_mesh.py index 1d000546..70bafcae 100644 --- a/lapy/tria_mesh.py +++ b/lapy/tria_mesh.py @@ -1850,6 +1850,152 @@ def __reduce_edges_to_path( else: return paths + def extract_level_paths( + self, + vfunc: np.ndarray, + level: float, + ) -> list[polygon.Polygon]: + """Extract level set paths as Polygon objects with triangle/edge metadata. + + For a given scalar function on mesh vertices, extracts all paths where + the function equals the specified level. Returns polygons with embedded + metadata about which triangles and edges were intersected. + Handles single open paths, closed loops, and multiple disconnected components. + + Parameters + ---------- + vfunc : np.ndarray + Scalar function values at vertices, shape (n_vertices,). Must be 1D. + level : float + Level set value to extract. + + Returns + ------- + list[polygon.Polygon] + List of Polygon objects, one for each connected level curve component. + Each polygon has the following additional attributes: + - points : np.ndarray of shape (n_points, 3) - 3D coordinates on level curve + (access via polygon.points or polygon.get_points()) + - closed : bool - whether the curve is closed + - tria_idx : np.ndarray of shape (n_segments,) - triangle index for each segment + - edge_vidx : np.ndarray of shape (n_points, 2) - mesh vertex indices for edge + - edge_bary : np.ndarray of shape (n_points,) - barycentric coordinate [0,1] + along edge where level set intersects (0=first vertex, 1=second vertex) + + Raises + ------ + ValueError + If vfunc is not 1-dimensional. + """ + if vfunc.ndim > 1: + raise ValueError(f"vfunc needs to be 1-dim, but is {vfunc.ndim}-dim!") + + # Get intersecting triangles + intersect = vfunc[self.t] > level + t_idx = np.where( + np.logical_or( + np.sum(intersect, axis=1) == 1, np.sum(intersect, axis=1) == 2 + ) + )[0] + + if t_idx.size == 0: + return [] + + # Reduce to triangles that intersect with level + t_level = self.t[t_idx, :] + intersect = intersect[t_idx, :] + + # Invert triangles with two true values so single vertex is true + intersect[np.sum(intersect, axis=1) > 1, :] = np.logical_not( + intersect[np.sum(intersect, axis=1) > 1, :] + ) + + # Get index within triangle with single vertex + idx_single = np.argmax(intersect, axis=1) + idx_o1 = (idx_single + 1) % 3 + idx_o2 = (idx_single + 2) % 3 + + # Get global indices + gidx0 = t_level[np.arange(t_level.shape[0]), idx_single] + gidx1 = t_level[np.arange(t_level.shape[0]), idx_o1] + gidx2 = t_level[np.arange(t_level.shape[0]), idx_o2] + + # Sort edge indices (rows are triangles, cols are the two vertex ids) + # This creates unique edge identifiers + gg1 = np.sort( + np.concatenate((gidx0[:, np.newaxis], gidx1[:, np.newaxis]), axis=1) + ) + gg2 = np.sort( + np.concatenate((gidx0[:, np.newaxis], gidx2[:, np.newaxis]), axis=1) + ) + + # Concatenate all edges and get unique ones + gg = np.concatenate((gg1, gg2), axis=0) + gg_unique = np.unique(gg, axis=0) + + # Generate level set intersection points for unique edges + # Barycentric coordinate (0=first vertex, 1=second vertex) + xl = ((level - vfunc[gg_unique[:, 0]]) + / (vfunc[gg_unique[:, 1]] - vfunc[gg_unique[:, 0]])) + + # 3D points on unique edges + p = ((1 - xl)[:, np.newaxis] * self.v[gg_unique[:, 0]] + + xl[:, np.newaxis] * self.v[gg_unique[:, 1]]) + + # Fill sparse matrix with new point indices (+1 to distinguish from zero) + a_mat = sparse.csc_matrix( + (np.arange(gg_unique.shape[0]) + 1, (gg_unique[:, 0], gg_unique[:, 1])) + ) + + # For each triangle, create edge pair via lookup in matrix + edge_idxs = (np.concatenate((a_mat[gg1[:, 0], gg1[:, 1]], + a_mat[gg2[:, 0], gg2[:, 1]]), axis=0).T - 1) + + # Reduce edges to paths (returns list of paths for multiple components) + paths, path_edge_idxs = self._TriaMesh__reduce_edges_to_path(edge_idxs, get_edge_idx=True) + + # Build polygon objects for each connected component + polygons = [] + for path_nodes, path_edge_idx in zip(paths, path_edge_idxs): + # Get 3D coordinates for this path + poly_v = p[path_nodes, :] + + # Remove duplicate vertices (when levelset hits a vertex almost perfectly) + dd = ((poly_v[0:-1, :] - poly_v[1:, :]) ** 2).sum(1) + dd = np.append(dd, 1) # Never delete last node + eps = 0.000001 + keep_ids = dd > eps + poly_v = poly_v[keep_ids, :] + + # Get triangle indices for each edge + poly_tria_idx = t_idx[path_edge_idx[keep_ids[:-1]]] + + # Get edge vertex indices + poly_edge_vidx = gg_unique[path_nodes[keep_ids], :] + + # Get barycentric coordinates + poly_edge_bary = xl[path_nodes[keep_ids]] + + # Create polygon with metadata + # Use closed=None to let Polygon auto-detect based on endpoints + # This will automatically remove duplicate endpoint if present + n_points_before = poly_v.shape[0] + poly = polygon.Polygon(poly_v, closed=None) + n_points_after = poly.points.shape[0] + + # If Polygon removed duplicate endpoint, adjust metadata + if n_points_after < n_points_before: + # Remove last entry from metadata to match polygon points + poly_edge_vidx = poly_edge_vidx[:n_points_after] + poly_edge_bary = poly_edge_bary[:n_points_after] + + poly.tria_idx = poly_tria_idx + poly.edge_vidx = poly_edge_vidx + poly.edge_bary = poly_edge_bary + + polygons.append(poly) + + return polygons def level_path( self, @@ -1867,10 +2013,23 @@ def level_path( The points are sorted and returned as an ordered array of 3D coordinates together with the length of the level set path. + This implementation uses extract_level_paths internally to compute the + level set, ensuring consistent handling of closed loops and metadata. + .. note:: - Only works for level sets that represent a single non-intersecting - path with exactly one start and one endpoint (e.g., not closed)! + Works for level sets that represent a single path or closed loop. + For closed loops, the first and last points are identical (duplicated) + so you can detect closure via ``np.allclose(path[0], path[-1])``. + For open paths, the path has two distinct endpoints. + This function is kept mainly for backward compatability. + + **For more advanced use cases, consider using extract_level_paths() directly:** + + - Multiple disconnected components: extract_level_pathss returns all components + - Custom resampling: Get Polygon objects and use Polygon.resample() directly + - Metadata analysis: Access triangle indices and edge information per component + - Closed loop detection: Polygon objects have a ``closed`` attribute Parameters ---------- @@ -1880,6 +2039,8 @@ def level_path( Level set value to extract. get_tria_idx : bool, default=False If True, also return array of triangle indices for each path segment. + For closed loops with n points (including duplicate), returns n-1 triangle + indices. For open paths with n points, returns n-1 triangle indices. get_edges : bool, default=False If True, also return arrays of vertex indices (i,j) and relative positions for each 3D point along the intersecting edge. @@ -1891,12 +2052,15 @@ def level_path( ------- path : np.ndarray Array of shape (n, 3) containing 3D coordinates of vertices on the - level path. + level path. For closed loops, the last point duplicates the first point, + so ``np.allclose(path[0], path[-1])`` will be True. length : float - Length of the level set. + Total length of the level set path. For closed loops, includes the + closing segment from last to first point. tria_idx : np.ndarray Array of triangle indices for each segment on the path, shape (n-1,) - if the path has length n. Only returned if get_tria_idx is True. + where n is the number of points (including duplicate for closed loops). + Only returned if get_tria_idx is True. edges_vidxs : np.ndarray Array of shape (n, 2) with vertex indices (i,j) for each 3D point, defining the vertices of the original mesh edge intersecting the @@ -1912,116 +2076,90 @@ def level_path( ------ ValueError If vfunc is not 1-dimensional. + If multiple disconnected level paths are found (use extract_level_paths). If n_points is combined with get_tria_idx=True. If n_points is combined with get_edges=True. + + See Also + -------- + extract_level_paths : Extract multiple disconnected level paths as Polygon objects. + Recommended for advanced use cases with multiple components, custom resampling, + or when you need explicit closed/open information. + + Examples + -------- + >>> # Extract a simple level curve + >>> path, length = mesh.level_path(vfunc, 0.5) + >>> print(f"Path has {path.shape[0]} points, length {length:.2f}") + + >>> # Check if the path is closed + >>> is_closed = np.allclose(path[0], path[-1]) + >>> print(f"Path is closed: {is_closed}") + + >>> # Get triangle indices for each segment + >>> path, length, tria_idx = mesh.level_path(vfunc, 0.5, get_tria_idx=True) + + >>> # Get complete edge metadata + >>> path, length, edges, relpos = mesh.level_path(vfunc, 0.5, get_edges=True) + + >>> # Resample to fixed number of points + >>> path, length = mesh.level_path(vfunc, 0.5, n_points=100) + + >>> # For multiple components, use extract_level_paths instead + >>> curves = mesh.extract_level_paths(vfunc, 0.5) + >>> for i, curve in enumerate(curves): + >>> print(f"Component {i}: {curve.points.shape[0]} points, closed={curve.closed}") """ - if vfunc.ndim > 1: - raise ValueError(f"vfunc needs to be 1-dim, but is {vfunc.ndim}-dim!") - # get intersecting triangles - intersect = vfunc[self.t] > level - t_idx = np.where( - np.logical_or( - np.sum(intersect, axis=1) == 1, np.sum(intersect, axis=1) == 2 - ) - )[0] - # reduce to triangles that intersect with level: - t_level = self.t[t_idx, :] - intersect = intersect[t_idx, :] - # trias have one vertex on one side and two on the other side of the level set - # invert trias with two true values, so that single vertex is true - intersect[np.sum(intersect, axis=1) > 1, :] = np.logical_not( - intersect[np.sum(intersect, axis=1) > 1, :] - ) - # get idx within tria with single vertex: - idx_single = np.argmax(intersect, axis=1) - idx_o1 = (idx_single + 1) % 3 - idx_o2 = (idx_single + 2) % 3 - # get global idx - gidx0 = t_level[np.arange(t_level.shape[0]), idx_single] - gidx1 = t_level[np.arange(t_level.shape[0]), idx_o1] - gidx2 = t_level[np.arange(t_level.shape[0]), idx_o2] - # sort edge indices (rows are trias, cols are the two vertex ids) - gg1 = np.sort( - np.concatenate((gidx0[:, np.newaxis], gidx1[:, np.newaxis]), axis=1) - ) - gg2 = np.sort( - np.concatenate((gidx0[:, np.newaxis], gidx2[:, np.newaxis]), axis=1) - ) - # concatenate all and get unique ones - gg = np.concatenate((gg1, gg2), axis=0) - gg_unique = np.unique(gg, axis=0) - # generate level set intersection points for unique edges - xl = ((level - vfunc[gg_unique[:, 0]]) - / ( vfunc[gg_unique[:, 1]] - vfunc[gg_unique[:, 0]])) - p = ((1 - xl)[:, np.newaxis] * self.v[gg_unique[:, 0]] - + xl[:, np.newaxis] * self.v[gg_unique[:, 1]]) - # fill sparse matrix with new point indices (+1 to distinguish from zero) - a_mat = sparse.csc_matrix( - (np.arange(gg_unique.shape[0]) + 1, (gg_unique[:, 0], gg_unique[:, 1])) - ) - # for each tria create one edge via lookup in matrix - edge_idxs = ( np.concatenate((a_mat[gg1[:, 0], gg1[:, 1]], - a_mat[gg2[:, 0], gg2[:, 1]]), axis=0).T - 1 ) - # lengths computation - p1 = np.squeeze(p[edge_idxs[:, 0]]) - p2 = np.squeeze(p[edge_idxs[:, 1]]) - llength = np.sqrt(((p1 - p2) ** 2).sum(1)).sum() - # compute path from unordered, not-directed edge list - # __reduce_edges_to_path now returns lists (can have multiple components) - if get_tria_idx: - paths, edge_idxs_list = TriaMesh.__reduce_edges_to_path( - edge_idxs, get_edge_idx=True + # Use extract_level_paths to get polygons + curves = self.extract_level_paths(vfunc, level) + + # level_path expects single component + if len(curves) != 1: + raise ValueError( + f"Found {len(curves)} disconnected level curves. " + f"Use extract_level_paths() for multiple components." ) - # level_path expects single component - if len(paths) != 1: - raise ValueError( - f"Found {len(paths)} disconnected level curves. " - f"Use extract_level_curves() for multiple components." - ) - path = paths[0] - edge_idx = edge_idxs_list[0] - # translate local edge id to global tria id - tria_idx = t_idx[edge_idx] - else: - paths = TriaMesh.__reduce_edges_to_path(edge_idxs, get_edge_idx=False) - # level_path expects single component - if len(paths) != 1: - raise ValueError( - f"Found {len(paths)} disconnected level curves. " - f"Use extract_level_curves() for multiple components." - ) - path = paths[0] - - # remove duplicate vertices (happens when levelset hits a vertex almost - # perfectly) - path3d = p[path, :] - dd = ((path3d[0:-1, :] - path3d[1:, :]) ** 2).sum(1) - # append 1 (never delete last node, if identical to the one before, we delete - # the one before) - dd = np.append(dd, 1) - eps = 0.000001 - keep_ids = dd > eps - path3d = path3d[keep_ids, :] - edges_vidxs = gg_unique[path, :] - edges_vidxs = edges_vidxs[keep_ids, :] - edges_relpos = xl[path] - edges_relpos = edges_relpos[keep_ids] - if get_tria_idx: - if n_points: + + # Get the single curve + curve = curves[0] + + # Extract data from polygon + path3d = curve.points + edges_vidxs = curve.edge_vidx + edges_relpos = curve.edge_bary + + # Compute length using polygon's length method + length = curve.length() + + # For closed loops, duplicate the last point to match first point + # This allows users to detect closure via np.allclose(path[0], path[-1]) + # and maintains backward compatibility with the original level_path behavior + if curve.closed: + path3d = np.vstack([path3d, path3d[0:1]]) + edges_vidxs = np.vstack([edges_vidxs, edges_vidxs[0:1]]) + edges_relpos = np.append(edges_relpos, edges_relpos[0]) + + # Handle optional resampling + if n_points: + if get_tria_idx: raise ValueError("n_points cannot be combined with get_tria_idx=True.") - tria_idx = tria_idx[dd[:-1] > eps] if get_edges: - return path3d, llength, tria_idx, edges_vidxs, edges_relpos + raise ValueError("n_points cannot be combined with get_edges=True.") + poly = polygon.Polygon(path3d, closed=curve.closed) + path3d = poly.resample(n_points=n_points, n_iter=3, inplace=False) + path3d = path3d.get_points() + + # Build return tuple based on options + if get_tria_idx: + tria_idx = curve.tria_idx + if get_edges: + return path3d, length, tria_idx, edges_vidxs, edges_relpos else: - return path3d, llength, tria_idx + return path3d, length, tria_idx else: - if n_points: - if get_edges: - raise ValueError("n_points cannot be combined with get_edges=True.") - poly = polygon.Polygon(path3d, closed=False) - path3d = poly.resample(n_points=n_points, n_iter=3, inplace=False) - path3d = path3d.get_points() if get_edges: - return path3d, llength, edges_vidxs, edges_relpos + return path3d, length, edges_vidxs, edges_relpos else: - return path3d, llength + return path3d, length + + diff --git a/lapy/utils/tests/test_tria_mesh.py b/lapy/utils/tests/test_tria_mesh.py index a3b0a4a5..7b3d57c1 100644 --- a/lapy/utils/tests/test_tria_mesh.py +++ b/lapy/utils/tests/test_tria_mesh.py @@ -630,3 +630,475 @@ def test_2d_mesh_support(): assert v2d.shape == (4, 2), "Should return 2D vertices" np.testing.assert_array_almost_equal(v2d, vertices_2d) + +def test_critical_points_simple(): + """ + Test critical_points on a simple mesh with known extrema. + """ + # Create a simple triangular mesh with 3 peaks and a valley + vertices = np.array([ + [0.0, 0.0, 0.5], # 0: center (middle height) + [1.0, 0.0, 1.0], # 1: peak + [0.0, 1.0, 0.0], # 2: valley + [-1.0, 0.0, 1.0], # 3: peak + ]) + triangles = np.array([ + [0, 1, 2], + [0, 2, 3], + [0, 3, 1], + ]) + mesh = TriaMesh(vertices, triangles) + + # Height function: z-coordinate + vfunc = vertices[:, 2] + + # Compute critical points + minima, maxima, saddles, saddle_orders = mesh.critical_points(vfunc) + + # Expected: vertex 2 is minimum (z=0), vertices 1 and 3 are maxima (z=1) + assert 2 in minima, f"Expected vertex 2 (valley) as minimum, got minima={minima}" + assert 1 in maxima or 3 in maxima, f"Expected vertices 1 or 3 (peaks) as maxima, got maxima={maxima}" + # Vertex 0 is at mid-height, surrounded by higher and lower neighbors - could be saddle or regular + + # Simpler test: just verify the function runs without error + assert len(minima) + len(maxima) + len(saddles) > 0, "Should find some critical points" + + +def test_critical_points_saddle(): + """ + Test critical_points on a mesh with a saddle point. + """ + # Create a simple saddle mesh: 3x3 grid in xy-plane + # Height creates saddle at center + vertices = np.array([ + [-1, -1, 1], # 0: corner (high) + [0, -1, 0], # 1: edge midpoint (low) + [1, -1, 1], # 2: corner (high) + [-1, 0, 0], # 3: edge midpoint (low) + [0, 0, 0.5], # 4: center (saddle) + [1, 0, 0], # 5: edge midpoint (low) + [-1, 1, 1], # 6: corner (high) + [0, 1, 0], # 7: edge midpoint (low) + [1, 1, 1], # 8: corner (high) + ], dtype=float) + + # Create triangular mesh from grid + triangles = np.array([ + [0, 1, 4], [0, 4, 3], # bottom-left quad + [1, 2, 5], [1, 5, 4], # bottom-right quad + [3, 4, 7], [3, 7, 6], # top-left quad + [4, 5, 8], [4, 8, 7], # top-right quad + ]) + mesh = TriaMesh(vertices, triangles) + + # Use z-coordinate as function + vfunc = vertices[:, 2] + + # Compute critical points + minima, maxima, saddles, saddle_orders = mesh.critical_points(vfunc) + + # Center should be a saddle (neighbors alternate high-low) + assert 4 in saddles, f"Expected vertex 4 (center) to be a saddle, saddles={saddles}" + # Check saddle order for center vertex + saddle_idx = np.where(saddles == 4)[0] + if len(saddle_idx) > 0: + order = saddle_orders[saddle_idx[0]] + assert order >= 2, f"Expected saddle order >= 2, got {order}" + + +def test_critical_points_with_ties(): + """ + Test critical_points with equal function values (tie-breaker rule). + """ + # Simple triangle mesh + vertices = np.array([ + [0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [0.5, 1.0, 0.0], + ]) + triangles = np.array([[0, 1, 2]]) + mesh = TriaMesh(vertices, triangles) + + # All vertices have same value (ties everywhere) + vfunc = np.array([1.0, 1.0, 1.0]) + + # Should not crash, tie-breaker should resolve ambiguities + minima, maxima, saddles, saddle_orders = mesh.critical_points(vfunc) + + # With tie-breaker, vertex with highest ID is treated as larger + # So vertex 2 should be maximum, vertex 0 should be minimum + assert 0 in minima, "Expected vertex 0 to be minimum with tie-breaker" + assert 2 in maxima, "Expected vertex 2 to be maximum with tie-breaker" + + +def test_critical_points_boundary(): + """ + Test critical_points on a mesh with boundary (non-closed). + """ + # Single triangle (has boundary) + vertices = np.array([ + [0.0, 0.0, 0.0], + [1.0, 0.0, 1.0], + [0.0, 1.0, 0.0], + ]) + triangles = np.array([[0, 1, 2]]) + mesh = TriaMesh(vertices, triangles) + + # Height function + vfunc = vertices[:, 2] + + minima, maxima, saddles, saddle_orders = mesh.critical_points(vfunc) + + # Vertex 1 is highest (z=1), vertices 0 and 2 are lowest (z=0) + assert 1 in maxima, f"Expected vertex 1 as maximum, got maxima={maxima}" + assert len(minima) >= 1, f"Expected at least one minimum, got {len(minima)}" + + +def test_extract_level_paths_simple(): + """ + Test extract_level_paths on a simple mesh with single level curve. + """ + # Create a simple mesh: unit square in xy-plane + vertices = np.array([ + [0.0, 0.0, 0.0], + [1.0, 0.0, 0.0], + [1.0, 1.0, 1.0], + [0.0, 1.0, 1.0], + ]) + triangles = np.array([ + [0, 1, 2], + [0, 2, 3], + ]) + mesh = TriaMesh(vertices, triangles) + + # Height function (z-coordinate) + vfunc = vertices[:, 2] + + # Extract level curve at z=0.5 + curves = mesh.extract_level_paths(vfunc, 0.5) + + # Should have at least one curve + assert len(curves) > 0, "Expected at least one level curve" + + # Check that returned objects are Polygon + from ...polygon import Polygon + for curve in curves: + assert isinstance(curve, Polygon), f"Expected Polygon, got {type(curve)}" + + # Check that polygons have required attributes + for curve in curves: + assert hasattr(curve, 'points'), "Polygon should have 'points' attribute" + assert hasattr(curve, 'closed'), "Polygon should have 'closed' attribute" + assert hasattr(curve, 'tria_idx'), "Polygon should have 'tria_idx' attribute" + assert hasattr(curve, 'edge_vidx'), "Polygon should have 'edge_vidx' attribute" + assert hasattr(curve, 'edge_bary'), "Polygon should have 'edge_bary' attribute" + + # Check that points are 3D + for curve in curves: + assert curve.points.shape[1] == 3, f"Expected 3D points, got shape {curve.points.shape}" + + +def test_extract_level_paths_closed_loop(): + """ + Test extract_level_paths on a mesh with closed loop level curve. + """ + # Create a pyramid mesh + vertices = np.array([ + [0.0, 0.0, 0.0], # base + [1.0, 0.0, 0.0], # base + [1.0, 1.0, 0.0], # base + [0.0, 1.0, 0.0], # base + [0.5, 0.5, 1.0], # apex + ]) + triangles = np.array([ + [0, 1, 4], + [1, 2, 4], + [2, 3, 4], + [3, 0, 4], + [0, 2, 1], + [0, 3, 2], + ]) + mesh = TriaMesh(vertices, triangles) + + # Height function + vfunc = vertices[:, 2] + + # Extract level curve at mid-height (should create closed loop around pyramid) + curves = mesh.extract_level_paths(vfunc, 0.5) + + assert len(curves) > 0, "Expected at least one level curve" + + # At least one curve should be closed + has_closed = any(curve.closed for curve in curves) + # Note: May or may not be closed depending on mesh topology + + # All points should be approximately at z=0.5 + for curve in curves: + z_coords = curve.points[:, 2] + np.testing.assert_allclose(z_coords, 0.5, atol=1e-5, + err_msg="Level curve points should be at z=0.5") + + +def test_extract_level_paths_no_intersection(): + """ + Test extract_level_paths when level doesn't intersect mesh. + """ + # Simple triangle + vertices = np.array([ + [0.0, 0.0, 0.0], + [1.0, 0.0, 1.0], + [0.0, 1.0, 1.0], + ]) + triangles = np.array([[0, 1, 2]]) + mesh = TriaMesh(vertices, triangles) + + # Height function + vfunc = vertices[:, 2] + + # Extract level curve at z=10.0 (way above mesh) + curves = mesh.extract_level_paths(vfunc, 10.0) + + # Should return empty list + assert len(curves) == 0, f"Expected no curves, got {len(curves)}" + + +def test_extract_level_paths_metadata(): + """ + Test that extract_level_paths returns correct metadata. + """ + # Create a simple mesh + vertices = np.array([ + [0.0, 0.0, 0.0], + [1.0, 0.0, 0.5], + [1.0, 1.0, 1.0], + [0.0, 1.0, 0.5], + ]) + triangles = np.array([ + [0, 1, 2], + [0, 2, 3], + ]) + mesh = TriaMesh(vertices, triangles) + + # Height function + vfunc = vertices[:, 2] + + # Extract level curve + curves = mesh.extract_level_paths(vfunc, 0.5) + + for curve in curves: + n_points = curve.points.shape[0] + + # Check tria_idx + assert hasattr(curve, 'tria_idx'), "Missing tria_idx" + # Number of segments (tria_idx) depends on whether curve is closed + if curve.closed: + expected_segments = n_points + else: + expected_segments = n_points - 1 + # Note: tria_idx may be shorter if duplicate points were removed + + # Check edge_vidx + assert hasattr(curve, 'edge_vidx'), "Missing edge_vidx" + assert curve.edge_vidx.shape == (n_points, 2), \ + f"edge_vidx should be (n_points, 2), got {curve.edge_vidx.shape}" + + # Check edge_bary + assert hasattr(curve, 'edge_bary'), "Missing edge_bary" + assert curve.edge_bary.shape == (n_points,), \ + f"edge_bary should be (n_points,), got {curve.edge_bary.shape}" + # Barycentric coordinates should be in [0, 1] + assert np.all(curve.edge_bary >= 0) and np.all(curve.edge_bary <= 1), \ + "edge_bary should be in [0, 1]" + + +def test_level_path(): + """ + Test level_path function for single path extraction with various options. + """ + # Create a simple mesh with a clear level curve + vertices = np.array([ + [0.0, 0.0, 0.0], + [1.0, 0.0, 0.5], + [1.0, 1.0, 1.0], + [0.0, 1.0, 0.5], + ]) + triangles = np.array([ + [0, 1, 2], + [0, 2, 3], + ]) + mesh = TriaMesh(vertices, triangles) + vfunc = vertices[:, 2] + level = 0.5 + + # Test 1: Basic usage - get path and length + path, length = mesh.level_path(vfunc, level) + + assert path.ndim == 2, "Path should be 2D array" + assert path.shape[1] == 3, f"Path should have 3D points, got shape {path.shape}" + assert path.shape[0] >= 2, "Path should have at least 2 points" + assert length > 0, f"Path length should be positive, got {length}" + + # Check all points are approximately at the level + z_coords = path[:, 2] + np.testing.assert_allclose(z_coords, level, atol=1e-5, + err_msg=f"All path points should be at z={level}") + + # Test 2: Check if path is closed (first == last point) + is_closed = np.allclose(path[0], path[-1]) + if is_closed: + print(f" Path is closed (first point == last point)") + else: + print(f" Path is open (endpoints differ)") + + # Test 3: Get triangle indices + path_with_tria, length_with_tria, tria_idx = mesh.level_path( + vfunc, level, get_tria_idx=True + ) + + np.testing.assert_array_equal(path_with_tria, path, + err_msg="Path should be same with get_tria_idx=True") + assert length_with_tria == length, "Length should be same" + assert tria_idx.ndim == 1, "tria_idx should be 1D array" + # For n points, we have n-1 segments + expected_tria_len = path.shape[0] - 1 + assert tria_idx.shape[0] == expected_tria_len, \ + f"tria_idx should have {expected_tria_len} elements, got {tria_idx.shape[0]}" + # Triangle indices should be valid + assert np.all(tria_idx >= 0), "Triangle indices should be non-negative" + assert np.all(tria_idx < len(triangles)), \ + f"Triangle indices should be < {len(triangles)}" + + # Test 4: Get edge information + path_with_edges, length_with_edges, edges_vidxs, edges_relpos = mesh.level_path( + vfunc, level, get_edges=True + ) + + np.testing.assert_array_equal(path_with_edges, path, + err_msg="Path should be same with get_edges=True") + assert length_with_edges == length, "Length should be same" + assert edges_vidxs.shape == (path.shape[0], 2), \ + f"edges_vidxs should be (n_points, 2), got {edges_vidxs.shape}" + assert edges_relpos.shape == (path.shape[0],), \ + f"edges_relpos should be (n_points,), got {edges_relpos.shape}" + # Barycentric coordinates should be in [0, 1] + assert np.all(edges_relpos >= 0) and np.all(edges_relpos <= 1), \ + "Barycentric coordinates should be in [0, 1]" + + # Test 5: Get both triangle and edge information + result = mesh.level_path(vfunc, level, get_tria_idx=True, get_edges=True) + path_full, length_full, tria_idx_full, edges_vidxs_full, edges_relpos_full = result + + np.testing.assert_array_equal(path_full, path, + err_msg="Path should be consistent") + np.testing.assert_array_equal(tria_idx_full, tria_idx, + err_msg="tria_idx should be consistent") + np.testing.assert_array_equal(edges_vidxs_full, edges_vidxs, + err_msg="edges_vidxs should be consistent") + np.testing.assert_array_equal(edges_relpos_full, edges_relpos, + err_msg="edges_relpos should be consistent") + + # Test 6: Resampling to fixed number of points + n_resample = 50 + path_resampled, length_resampled = mesh.level_path( + vfunc, level, n_points=n_resample + ) + + assert path_resampled.shape[0] == n_resample, \ + f"Resampled path should have {n_resample} points, got {path_resampled.shape[0]}" + assert path_resampled.shape[1] == 3, "Resampled path should have 3D points" + # Length should be approximately the same + np.testing.assert_allclose(length_resampled, length, rtol=0.1, + err_msg="Resampled length should be close to original") + + # Test 7: Error case - resampling with get_tria_idx should raise ValueError + try: + mesh.level_path(vfunc, level, get_tria_idx=True, n_points=50) + assert False, "Should raise ValueError when combining n_points with get_tria_idx" + except ValueError as e: + assert "n_points cannot be combined with get_tria_idx" in str(e) + + # Test 8: Error case - resampling with get_edges should raise ValueError + try: + mesh.level_path(vfunc, level, get_edges=True, n_points=50) + assert False, "Should raise ValueError when combining n_points with get_edges" + except ValueError as e: + assert "n_points cannot be combined with get_edges" in str(e) + + +def test_level_path_closed_loop(): + """ + Test level_path on a mesh that produces a closed loop. + """ + # Create a pyramid mesh + vertices = np.array([ + [0.0, 0.0, 0.0], # base + [1.0, 0.0, 0.0], # base + [1.0, 1.0, 0.0], # base + [0.0, 1.0, 0.0], # base + [0.5, 0.5, 1.0], # apex + ]) + triangles = np.array([ + [0, 1, 4], + [1, 2, 4], + [2, 3, 4], + [3, 0, 4], + [0, 2, 1], + [0, 3, 2], + ]) + mesh = TriaMesh(vertices, triangles) + vfunc = vertices[:, 2] + + # Extract level at mid-height (should create closed loop around pyramid) + path, length = mesh.level_path(vfunc, 0.5) + + # For closed loops, first and last points should be identical + is_closed = np.allclose(path[0], path[-1]) + assert is_closed, "Closed loop should have first point equal to last point" + + # All points should be at z=0.5 + np.testing.assert_allclose(path[:, 2], 0.5, atol=1e-5, + err_msg="All points should be at z=0.5") + + # Length should be positive + assert length > 0, "Closed loop should have positive length" + + +def test_level_path_multiple_components_error(): + """ + Test that level_path raises ValueError when multiple components are found. + """ + # Create a mesh that could have multiple level curves + # Two separate triangles at different locations + vertices = np.array([ + # First triangle + [0.0, 0.0, 0.0], + [1.0, 0.0, 1.0], + [0.5, 1.0, 0.5], + # Second triangle (disconnected) + [3.0, 0.0, 0.0], + [4.0, 0.0, 1.0], + [3.5, 1.0, 0.5], + ]) + triangles = np.array([ + [0, 1, 2], + [3, 4, 5], + ]) + mesh = TriaMesh(vertices, triangles) + vfunc = vertices[:, 2] + + # Try to extract level that intersects both triangles + # This should raise ValueError because level_path expects single component + try: + path, length = mesh.level_path(vfunc, 0.5) + # If we get here, either there's only 1 component or the test needs adjustment + # Check via extract_level_paths + curves = mesh.extract_level_paths(vfunc, 0.5) + if len(curves) > 1: + assert False, "Should have raised ValueError for multiple components" + # If only 1 component, that's fine - test is conditional + except ValueError as e: + # Expected error + assert "disconnected" in str(e).lower() or "multiple" in str(e).lower(), \ + f"Error message should mention multiple components: {e}" + + From 74b53dd1780d46668c3ccd5d1cf679710765b7c4 Mon Sep 17 00:00:00 2001 From: Martin Reuter Date: Wed, 21 Jan 2026 18:06:59 +0100 Subject: [PATCH 4/8] fix tests --- lapy/utils/tests/test_tria_mesh.py | 35 +++++++++--------------------- 1 file changed, 10 insertions(+), 25 deletions(-) diff --git a/lapy/utils/tests/test_tria_mesh.py b/lapy/utils/tests/test_tria_mesh.py index 7b3d57c1..557ec11c 100644 --- a/lapy/utils/tests/test_tria_mesh.py +++ b/lapy/utils/tests/test_tria_mesh.py @@ -826,10 +826,11 @@ def test_extract_level_paths_closed_loop(): # Extract level curve at mid-height (should create closed loop around pyramid) curves = mesh.extract_level_paths(vfunc, 0.5) - assert len(curves) > 0, "Expected at least one level curve" + assert len(curves) == 1, "Expected one level curve" + + # Curve should be closed + assert curves[0].closed, "Expected closed level curve" - # At least one curve should be closed - has_closed = any(curve.closed for curve in curves) # Note: May or may not be closed depending on mesh topology # All points should be approximately at z=0.5 @@ -891,11 +892,9 @@ def test_extract_level_paths_metadata(): # Check tria_idx assert hasattr(curve, 'tria_idx'), "Missing tria_idx" # Number of segments (tria_idx) depends on whether curve is closed - if curve.closed: - expected_segments = n_points - else: - expected_segments = n_points - 1 - # Note: tria_idx may be shorter if duplicate points were removed + expected_tria_len = n_points if curve.closed else n_points - 1 + assert curve.tria_idx.shape == (expected_tria_len,), \ + f"tria_idx should be ({expected_tria_len},), got {curve.tria_idx.shape}" # Check edge_vidx assert hasattr(curve, 'edge_vidx'), "Missing edge_vidx" @@ -946,9 +945,9 @@ def test_level_path(): # Test 2: Check if path is closed (first == last point) is_closed = np.allclose(path[0], path[-1]) if is_closed: - print(f" Path is closed (first point == last point)") + print(" Path is closed (first point == last point)") else: - print(f" Path is open (endpoints differ)") + print(" Path is open (endpoints differ)") # Test 3: Get triangle indices path_with_tria, length_with_tria, tria_idx = mesh.level_path( @@ -1010,20 +1009,6 @@ def test_level_path(): np.testing.assert_allclose(length_resampled, length, rtol=0.1, err_msg="Resampled length should be close to original") - # Test 7: Error case - resampling with get_tria_idx should raise ValueError - try: - mesh.level_path(vfunc, level, get_tria_idx=True, n_points=50) - assert False, "Should raise ValueError when combining n_points with get_tria_idx" - except ValueError as e: - assert "n_points cannot be combined with get_tria_idx" in str(e) - - # Test 8: Error case - resampling with get_edges should raise ValueError - try: - mesh.level_path(vfunc, level, get_edges=True, n_points=50) - assert False, "Should raise ValueError when combining n_points with get_edges" - except ValueError as e: - assert "n_points cannot be combined with get_edges" in str(e) - def test_level_path_closed_loop(): """ @@ -1094,7 +1079,7 @@ def test_level_path_multiple_components_error(): # Check via extract_level_paths curves = mesh.extract_level_paths(vfunc, 0.5) if len(curves) > 1: - assert False, "Should have raised ValueError for multiple components" + assert len(curves) == 1, "Should have raised ValueError for multiple components" # If only 1 component, that's fine - test is conditional except ValueError as e: # Expected error From 2c00321c7a116eb977e3f66514727c232171cca5 Mon Sep 17 00:00:00 2001 From: Martin Reuter Date: Wed, 21 Jan 2026 18:24:25 +0100 Subject: [PATCH 5/8] fix spacing etc --- lapy/tria_mesh.py | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/lapy/tria_mesh.py b/lapy/tria_mesh.py index 70bafcae..27eeb3ec 100644 --- a/lapy/tria_mesh.py +++ b/lapy/tria_mesh.py @@ -1065,13 +1065,11 @@ def critical_points(self, vfunc: np.ndarray) -> tuple[np.ndarray, np.ndarray, np 2. Minima: all differences are positive (all neighbors higher) 3. Maxima: all differences are negative (all neighbors lower) 4. Saddles: count sign flips across opposite edges in triangles at vertex - - Regular point: 2 sign flips - - Simple saddle (order 2): 4 sign flips - - Higher-order saddle (order n): 2n sign flips, order = n_flips / 2 + - Regular point: 2 sign flips + - Simple saddle (order 2): 4 sign flips + - Higher-order saddle (order n): 2n sign flips, order = n_flips / 2 5. Tie-breaker: when two vertices have equal function values, the vertex - with the higher vertex ID is treated as slightly larger to remove ambiguity. - - The implementation is fully vectorized using sparse matrices and numpy operations. + with the higher vertex ID is treated as slightly larger to remove ambiguity. """ vfunc = np.asarray(vfunc) if len(vfunc) != self.v.shape[0]: @@ -1875,12 +1873,11 @@ def extract_level_paths( List of Polygon objects, one for each connected level curve component. Each polygon has the following additional attributes: - points : np.ndarray of shape (n_points, 3) - 3D coordinates on level curve - (access via polygon.points or polygon.get_points()) - closed : bool - whether the curve is closed - tria_idx : np.ndarray of shape (n_segments,) - triangle index for each segment - edge_vidx : np.ndarray of shape (n_points, 2) - mesh vertex indices for edge - edge_bary : np.ndarray of shape (n_points,) - barycentric coordinate [0,1] - along edge where level set intersects (0=first vertex, 1=second vertex) + along edge where level set intersects (0=first vertex, 1=second vertex) Raises ------ @@ -2022,7 +2019,7 @@ def level_path( For closed loops, the first and last points are identical (duplicated) so you can detect closure via ``np.allclose(path[0], path[-1])``. For open paths, the path has two distinct endpoints. - This function is kept mainly for backward compatability. + This function is kept mainly for backward compatibility. **For more advanced use cases, consider using extract_level_paths() directly:** From cc9a49cb0afc74e9f6ddeb78d095ab36a3340abe Mon Sep 17 00:00:00 2001 From: Martin Reuter Date: Wed, 21 Jan 2026 18:25:21 +0100 Subject: [PATCH 6/8] Update lapy/tria_mesh.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- lapy/tria_mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lapy/tria_mesh.py b/lapy/tria_mesh.py index 27eeb3ec..4be16d87 100644 --- a/lapy/tria_mesh.py +++ b/lapy/tria_mesh.py @@ -2023,7 +2023,7 @@ def level_path( **For more advanced use cases, consider using extract_level_paths() directly:** - - Multiple disconnected components: extract_level_pathss returns all components + - Multiple disconnected components: extract_level_paths returns all components - Custom resampling: Get Polygon objects and use Polygon.resample() directly - Metadata analysis: Access triangle indices and edge information per component - Closed loop detection: Polygon objects have a ``closed`` attribute From e05d4101c443105298f5ae3d746247923c0d08e7 Mon Sep 17 00:00:00 2001 From: Martin Reuter Date: Wed, 21 Jan 2026 18:41:30 +0100 Subject: [PATCH 7/8] fix spacing etc --- lapy/tria_mesh.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/lapy/tria_mesh.py b/lapy/tria_mesh.py index 4be16d87..43d34549 100644 --- a/lapy/tria_mesh.py +++ b/lapy/tria_mesh.py @@ -1061,15 +1061,18 @@ def critical_points(self, vfunc: np.ndarray) -> tuple[np.ndarray, np.ndarray, np Notes ----- The algorithm works by: + 1. For each vertex, compute difference: neighbor_value - vertex_value 2. Minima: all differences are positive (all neighbors higher) 3. Maxima: all differences are negative (all neighbors lower) 4. Saddles: count sign flips across opposite edges in triangles at vertex - - Regular point: 2 sign flips - - Simple saddle (order 2): 4 sign flips - - Higher-order saddle (order n): 2n sign flips, order = n_flips / 2 + + - Regular point: 2 sign flips + - Simple saddle (order 2): 4 sign flips + - Higher-order saddle (order n): 2n sign flips, order = n_flips / 2 + 5. Tie-breaker: when two vertices have equal function values, the vertex - with the higher vertex ID is treated as slightly larger to remove ambiguity. + with the higher vertex ID is treated as slightly larger to remove ambiguity. """ vfunc = np.asarray(vfunc) if len(vfunc) != self.v.shape[0]: @@ -1872,12 +1875,13 @@ def extract_level_paths( list[polygon.Polygon] List of Polygon objects, one for each connected level curve component. Each polygon has the following additional attributes: + - points : np.ndarray of shape (n_points, 3) - 3D coordinates on level curve - closed : bool - whether the curve is closed - tria_idx : np.ndarray of shape (n_segments,) - triangle index for each segment - edge_vidx : np.ndarray of shape (n_points, 2) - mesh vertex indices for edge - edge_bary : np.ndarray of shape (n_points,) - barycentric coordinate [0,1] - along edge where level set intersects (0=first vertex, 1=second vertex) + along edge where level set intersects (0=first vertex, 1=second vertex) Raises ------ From 0e39040b58da04aafb71220e5f98162a0fb281b9 Mon Sep 17 00:00:00 2001 From: Martin Reuter Date: Wed, 21 Jan 2026 18:57:24 +0100 Subject: [PATCH 8/8] remove some tests and fix loop closing after resample --- lapy/tria_mesh.py | 2 ++ lapy/utils/tests/test_tria_mesh.py | 40 ------------------------------ 2 files changed, 2 insertions(+), 40 deletions(-) diff --git a/lapy/tria_mesh.py b/lapy/tria_mesh.py index 43d34549..0da29ae5 100644 --- a/lapy/tria_mesh.py +++ b/lapy/tria_mesh.py @@ -2149,6 +2149,8 @@ def level_path( poly = polygon.Polygon(path3d, closed=curve.closed) path3d = poly.resample(n_points=n_points, n_iter=3, inplace=False) path3d = path3d.get_points() + if curve.closed: + path3d = np.vstack([path3d, path3d[0:1]]) # Build return tuple based on options if get_tria_idx: diff --git a/lapy/utils/tests/test_tria_mesh.py b/lapy/utils/tests/test_tria_mesh.py index 557ec11c..3811caf8 100644 --- a/lapy/utils/tests/test_tria_mesh.py +++ b/lapy/utils/tests/test_tria_mesh.py @@ -1047,43 +1047,3 @@ def test_level_path_closed_loop(): # Length should be positive assert length > 0, "Closed loop should have positive length" - -def test_level_path_multiple_components_error(): - """ - Test that level_path raises ValueError when multiple components are found. - """ - # Create a mesh that could have multiple level curves - # Two separate triangles at different locations - vertices = np.array([ - # First triangle - [0.0, 0.0, 0.0], - [1.0, 0.0, 1.0], - [0.5, 1.0, 0.5], - # Second triangle (disconnected) - [3.0, 0.0, 0.0], - [4.0, 0.0, 1.0], - [3.5, 1.0, 0.5], - ]) - triangles = np.array([ - [0, 1, 2], - [3, 4, 5], - ]) - mesh = TriaMesh(vertices, triangles) - vfunc = vertices[:, 2] - - # Try to extract level that intersects both triangles - # This should raise ValueError because level_path expects single component - try: - path, length = mesh.level_path(vfunc, 0.5) - # If we get here, either there's only 1 component or the test needs adjustment - # Check via extract_level_paths - curves = mesh.extract_level_paths(vfunc, 0.5) - if len(curves) > 1: - assert len(curves) == 1, "Should have raised ValueError for multiple components" - # If only 1 component, that's fine - test is conditional - except ValueError as e: - # Expected error - assert "disconnected" in str(e).lower() or "multiple" in str(e).lower(), \ - f"Error message should mention multiple components: {e}" - -