lapy.TriaMesh¶

class lapy.TriaMesh(v, t, fsinfo=None)[source]¶

Class representing a triangle mesh.

This is an efficient implementation of a triangle mesh data structure with core functionality using sparse matrices internally (Scipy).

Parameters:

varray_like: List of lists of 2 or 3 float coordinates. For 2D vertices (n, 2), they will be automatically padded with z=0 to make them 3D internally. Must not be empty.
tarray_like: List of lists of 3 int of indices (>= 0) into v array Ordering is important: All triangles should be oriented in the same way (counter-clockwise, when looking from above). Must not be empty.
fsinfodict | None: FreeSurfer Surface Header Info.

Attributes

v	(array_like) List of lists of 3 float coordinates (internally always 3D).
t	(array_like) List of lists of 3 int of indices (>= 0) into `v` array.
adj_sym	(csc_matrix) Symmetric adjacency matrix as csc sparse matrix.
adj_dir	(csc_matrix) Directed adjacency matrix as csc sparse matrix.
fsinfo	(dict \| None) FreeSurfer Surface Header Info.
_is_2d	(bool) Internal flag indicating if the mesh was created with 2D vertices.

Methods

`area`()	Compute the total surface area of triangle mesh.
`avg_edge_length`()	Compute the average edge length of the mesh.
`boundary_loops`()	Compute boundary loops of the mesh.
`centroid`()	Compute centroid of triangle mesh as a weighted average of triangle centers.
`connected_components`()	Compute connected components of the mesh.
`construct_adj_dir_tidx`()	Construct directed adjacency matrix (edge graph) of triangle mesh.
`critical_points`(vfunc)	Compute critical points (extrema and saddles) of a function on mesh vertices.
`curvature`([smoothit])	Compute various curvature values at vertices.
`curvature_tria`([smoothit])	Compute min and max curvature and directions on triangle faces.
`edges`([with_boundary])	Compute vertices and adjacent triangle ids for each edge.
`euler`()	Compute the Euler Characteristic.
`extract_level_paths`(vfunc, level[, eps])	Extract level set paths as Polygon objects with triangle/edge metadata.
`get_vertices`([original_dim])	Get mesh vertices.
`has_free_vertices`()	Check if the vertex list has more vertices than what is used in triangles.
`is_2d`()	Check if the mesh was created with 2D vertices.
`is_boundary`()	Check which vertices are on the boundary.
`is_closed`()	Check if triangle mesh is closed (no boundary edges).
`is_manifold`()	Check if triangle mesh is manifold (no edges with >2 triangles).
`is_oriented`()	Check if triangle mesh is oriented.
`keep_largest_connected_component_`([clean])	Keep only the largest connected component of the mesh.
`level_length`(vfunc, level)	Compute the length of level sets.
`level_path`(vfunc, level[, get_tria_idx, ...])	Extract a single level set path as an array of 3D points.
`map_tfunc_to_vfunc`(tfunc[, weighted])	Map triangle function to vertices by attributing 1/3 to each vertex.
`map_vfunc_to_tfunc`(vfunc)	Map vertex function to triangles by averaging vertex values.
`normal_offset_`(d)	Move vertices along their normals by distance d.
`normalize_`()	Normalize TriaMesh to unit surface area and centroid at the origin.
`orient_`()	Re-orient triangles of manifold mesh to be consistent.
`read_fssurf`(filename)	Load triangle mesh from FreeSurfer surface geometry file.
`read_gifti`(filename)	Load triangle mesh from a GIFTI surface file.
`read_off`(filename)	Load triangle mesh from OFF txt file.
`read_vtk`(filename)	Load triangle mesh from VTK txt file.
`refine_`([it])	Refine the triangle mesh by placing new vertex on each edge midpoint.
`rm_free_vertices_`()	Remove unused (free) vertices.
`smooth_`([n])	Smooth the mesh iteratively in-place using Taubin smoothing.
`smooth_laplace`([vfunc, n, lambda_, mat])	Smooth the mesh or a vertex function using Laplace smoothing.
`smooth_taubin`([vfunc, n, lambda_, mu])	Smooth the mesh or a vertex function using Taubin smoothing.
`smooth_vfunc`(vfunc[, n])	Smooth the mesh or a vertex function iteratively.
`tria_areas`()	Compute the area of triangles using Heron's formula.
`tria_normals`()	Compute triangle normals.
`tria_qualities`()	Compute triangle quality for each triangle in mesh.
`vertex_areas`()	Compute the area associated to each vertex (1/3 of one-ring trias).
`vertex_degrees`()	Compute the vertex degrees (number of edges at each vertex).
`vertex_normals`()	Compute vertex normals.
`volume`()	Compute the volume of closed triangle mesh, summing tetrahedra at origin.
`write_fssurf`(filename[, image, ...])	Save as Freesurfer Surface Geometry file (wrap Nibabel).
`write_gifti`(filename)	Save mesh as a GIFTI surface file.
`write_vtk`(filename)	Save mesh as VTK file.

Raises:

ValueError: If mesh has no triangles or vertices (empty mesh). If vertices don’t have 2 or 3 coordinates. If triangles don’t have exactly 3 vertices. If triangle indices exceed number of vertices.

Notes

The class has static class methods to read triangle meshes from FreeSurfer, OFF, and VTK file formats.

When 2D vertices are provided, they are internally padded with z=0 to create 3D vertices. This allows all geometric operations to work correctly while maintaining compatibility with 2D mesh data.

Empty meshes are not allowed. Use class methods (read_fssurf, read_vtk, read_off, read_gifti) to load mesh data from files.

classmethod read_fssurf(filename)[source]¶

Load triangle mesh from FreeSurfer surface geometry file.

Parameters:

filenamestr: Filename to load, supporting .pial, .white, .sphere etc.

Returns:

TriaMesh: Loaded triangle mesh.

classmethod read_gifti(filename)[source]¶

Load triangle mesh from a GIFTI surface file.

GIFTI (.surf.gii) is the standard surface format used by HCP, FSL, FreeSurfer (since v6), and many other neuroimaging tools.

Parameters:

filenamestr: Path to a GIFTI surface file (e.g. lh.pial.surf.gii).

Returns:

TriaMesh: Loaded triangle mesh. Vertex coordinates are in the coordinate system stored in the file (usually surface RAS or MNI).

Raises:

ImportError: If nibabel is not installed.
OSError: If the file is not found or not readable.
ValueError: If the GIFTI file does not contain both a POINTSET and a TRIANGLE data array.

Parameters:

filename (str)

Return type:

TriaMesh

Notes

The GIFTI format supports multiple coordinate systems per file. This function uses whatever coordinate system nibabel exposes by default, which corresponds to the first coordinate system listed in the file (typically surface RAS for FreeSurfer output).

Examples

>>> from lapy import TriaMesh
>>> tria = TriaMesh.read_gifti("lh.pial.surf.gii")

classmethod read_off(filename)[source]¶

Load triangle mesh from OFF txt file.

Parameters:

filenamestr: Filename to load.

Returns:

TriaMesh: Loaded triangle mesh.

classmethod read_vtk(filename)[source]¶

Load triangle mesh from VTK txt file.

Parameters:

filenamestr: Filename to load.

Returns:

TriaMesh: Loaded triangle mesh.

area()[source]¶

Compute the total surface area of triangle mesh.

Returns:

float: Total surface area (sum of all triangle areas).

Return type:

float

avg_edge_length()[source]¶

Compute the average edge length of the mesh.

Returns:

float: Average edge length.

boundary_loops()[source]¶

Compute boundary loops of the mesh.

Meshes can have 0 or more boundary loops, which are cycles in the directed adjacency graph of the boundary edges. The mesh must be manifold and oriented.

Note

Could fail if loops are connected via a single vertex (like a figure 8). That case needs debugging.

Returns:

list of list: List of lists with boundary loops. Each inner list contains vertex indices forming a closed loop. Empty list if mesh is closed.

Raises:

ValueError: If mesh is not manifold (edges with more than 2 triangles). If mesh is not oriented.

Return type:

list

centroid()[source]¶

Compute centroid of triangle mesh as a weighted average of triangle centers.

The weight is determined by the triangle area. (This could be done much faster if a FEM lumped mass matrix M is already available where this would be M*v, because it is equivalent with averaging vertices weighted by vertex area)

Returns:

centroidnp.ndarray: The centroid coordinates of the mesh, shape (3,).
totalareafloat: The total area of the mesh.

Return type:

tuple[ndarray, float]

connected_components()[source]¶

Compute connected components of the mesh.

Uses scipy’s connected_components on the symmetric adjacency matrix.

Returns:

n_componentsint: Number of connected components.
labelsnp.ndarray: Label array of shape (n_vertices,) where labels[i] is the component ID of vertex i. Component IDs are integers from 0 to n_components-1.

Return type:

tuple[int, ndarray]

construct_adj_dir_tidx()[source]¶

Construct directed adjacency matrix (edge graph) of triangle mesh.

The directed adjacency matrix contains the triangle indices (only for non-manifold meshes). Operates only on triangles.

Returns:

csc_matrix: Similar to adj_dir, but stores the tria idx+1 instead of one in the matrix (allows lookup of vertex to tria).

critical_points(vfunc)[source]¶

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:

vfuncnp.ndarray: Real-valued function defined on mesh vertices, shape (n_vertices,).

Returns:

minimanp.ndarray: Array of vertex indices that are local minima, shape (n_minima,).
maximanp.ndarray: Array of vertex indices that are local maxima, shape (n_maxima,).
saddlesnp.ndarray: Array of vertex indices that are saddles (all orders), shape (n_saddles,).
saddle_ordersnp.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. Computed as ceil(n_flips / 2); for boundary saddles with an odd flip count the ceiling accounts for the implicit flip across the boundary.

Parameters:

vfunc (ndarray)

Return type:

tuple[ndarray, ndarray, ndarray, ndarray]

Notes

The algorithm works by:

For each vertex, compute difference: neighbor_value - vertex_value
Minima: all differences are positive (all neighbors higher)
Maxima: all differences are negative (all neighbors lower)
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
Tie-breaker: when two vertices have equal function values, the vertex with the higher vertex ID is treated as slightly larger to remove ambiguity.

curvature(smoothit=3)[source]¶

Compute various curvature values at vertices.

Computes principal curvature directions and values, mean and Gaussian curvature at each vertex using the anisotropic polygonal remeshing approach.

Note

For the algorithm see: Pierre Alliez, David Cohen-Steiner, Olivier Devillers, Bruno Levy, and Mathieu Desbrun. Anisotropic Polygonal Remeshing. ACM Transactions on Graphics, 2003.

Parameters:

smoothitint, default=3: Number of smoothing iterations on vertex functions.

Returns:

u_minnp.ndarray: Minimal curvature directions, shape (n_vertices, 3).
u_maxnp.ndarray: Maximal curvature directions, shape (n_vertices, 3).
c_minnp.ndarray: Minimal curvature values, shape (n_vertices,).
c_maxnp.ndarray: Maximal curvature values, shape (n_vertices,).
c_meannp.ndarray: Mean curvature (c_min + c_max) / 2.0, shape (n_vertices,).
c_gaussnp.ndarray: Gauss curvature c_min * c_max, shape (n_vertices,).
normalsnp.ndarray: Vertex normals, shape (n_vertices, 3).

Parameters:

smoothit (int)

Return type:

tuple[ndarray, ndarray, ndarray, ndarray, ndarray, ndarray, ndarray]

curvature_tria(smoothit=3)[source]¶

Compute min and max curvature and directions on triangle faces.

First computes curvature values on vertices and smooths them. Then maps them to triangles (averaging) and projects onto the triangle plane, ensuring orthogonality.

Parameters:

smoothitint, default=3: Number of smoothing iterations for curvature computation on vertices.

Returns:

u_minnp.ndarray: Minimal curvature direction on triangles, shape (n_triangles, 3).
u_maxnp.ndarray: Maximal curvature direction on triangles, shape (n_triangles, 3). Orthogonal to u_min and triangle normal.
c_minnp.ndarray: Minimal curvature values on triangles, shape (n_triangles,).
c_maxnp.ndarray: Maximal curvature values on triangles, shape (n_triangles,).

Parameters:

smoothit (int)

Return type:

tuple[ndarray, ndarray, ndarray, ndarray]

edges(with_boundary=False)[source]¶

Compute vertices and adjacent triangle ids for each edge.

Parameters:

with_boundarybool, default=False: If True, also return boundary half edges. If False, only return interior edges.

Returns:

vidsnp.ndarray: Array of shape (n_edges, 2) with starting and end vertex for each unique inner edge.
tidsnp.ndarray: Array of shape (n_edges, 2) with triangle containing the half edge from vids[:,0] to vids[:,1] in first column and the neighboring triangle in the second column.
bdrvidsnp.ndarray: If with_boundary is True: Array of shape (n_boundary_edges, 2) with each boundary half-edge. Only returned if with_boundary=True.
bdrtidsnp.ndarray: If with_boundary is True: Array of shape (n_boundary_edges, 1) with the associated triangle to each boundary edge. Only returned if with_boundary=True.

Raises:

ValueError: If mesh is not oriented.

Parameters:

with_boundary (bool)

Return type:

tuple[ndarray, …]

euler()[source]¶

Compute the Euler Characteristic.

The Euler characteristic is the number of vertices minus the number of edges plus the number of triangles (= #V - #E + #T). For example, it is 2 for the sphere and 0 for the torus. This operates only on triangles array.

Returns:

int: Euler characteristic.

Return type:

int

extract_level_paths(vfunc, level, eps=None)[source]¶

Extract level set paths as Polygon objects with triangle/edge metadata.

Extracts all paths where vfunc equals level on the mesh surface. Handles open paths, closed loops, multiple disconnected components, and level sets that pass exactly through mesh vertices:

Regular vertex: the two adjacent fragments are reconnected through the vertex into a single continuous path or loop.
Saddle vertex: the level set splits; each branch is emitted as an independent Polygon (potentially closed) and an info message is logged.
Boundary vertex: the open fragment is extended to include the vertex as its endpoint.

Parameters:

vfuncnp.ndarray: Scalar function values at vertices, shape (n_vertices,). Must be 1D.
levelfloat: Level set value to extract.
epsfloat or None, default=None: Tolerance for detecting vertices exactly on the level set. If None, defaults to max(1e-10, 1e-7 * (|level| + 1)).

Returns:

list[polygon.Polygon]

One Polygon per connected level curve component. Each polygon carries extra attributes set after construction:

tria_idx : np.ndarray, shape (n_segments,) — triangle index per segment. For closed curves n_segments == n_points; for open curves n_segments == n_points - 1.
edge_vidx : np.ndarray, shape (n_points, 2) — mesh vertex indices of the intersected edge per point. For special-vertex points both indices are equal (degenerate edge).
edge_bary : np.ndarray, shape (n_points,) — barycentric coordinate in [0, 1] along each edge (0 = first vertex, 1 = second vertex); 0.0 for special-vertex points.

Raises:

ValueError: If vfunc is not 1-dimensional.

Parameters:

vfunc (ndarray)
level (float)
eps (float | None)

Return type:

list[Polygon]

Notes

When two or more vertices of the same triangle lie exactly on the level set (level set runs along a mesh edge), a warning is issued and those triangles are skipped. The returned curves may be incomplete.

get_vertices(original_dim=False)[source]¶

Get mesh vertices.

Parameters:

original_dimbool, default=False: If True and mesh was created with 2D vertices, return vertices in original 2D format (without z-coordinate). If False, always return 3D vertices.

Returns:

np.ndarray: Vertex array of shape (n, 2) or (n, 3) depending on original_dim.

Parameters:

original_dim (bool)

Return type:

ndarray

has_free_vertices()[source]¶

Check if the vertex list has more vertices than what is used in triangles.

Free vertices are vertices that are not referenced by any triangle.

Returns:

bool: True if vertex list has more vertices than used in triangles, False otherwise.

Return type:

bool

is_2d()[source]¶

Check if the mesh was created with 2D vertices.

Returns:

bool: True if mesh was created with 2D vertices, False otherwise.

Return type:

bool

is_boundary()[source]¶

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).

Return type:

ndarray

is_closed()[source]¶

Check if triangle mesh is closed (no boundary edges).

Operates only on triangles

Returns:

bool: True if no boundary edges in adj matrix.

is_manifold()[source]¶

Check if triangle mesh is manifold (no edges with >2 triangles).

Operates only on triangles

Returns:

bool: True if no edges with > 2 triangles.

is_oriented()[source]¶

Check if triangle mesh is oriented.

True if all triangles are oriented counter-clockwise, when looking from above. Operates only on triangles.

Returns:

bool: True if max(adj_directed)=1.

Return type:

bool

keep_largest_connected_component_(clean=True)[source]¶

Keep only the largest connected component of the mesh.

Modifies the mesh in-place.

Parameters:

cleanbool, default=True: If True, remove vertices that are not used in the largest component. If False, vertices are kept but triangles from other components are removed, creating free vertices.

Returns:

vkeepnp.ndarray or None: Indices (from original list) of kept vertices if clean=True, else None.
vdelnp.ndarray or None: Indices of deleted (unused) vertices if clean=True, else None.

Parameters:

clean (bool)

Return type:

tuple[ndarray | None, ndarray | None]

level_length(vfunc, level)[source]¶

Compute the length of level sets.

For a scalar function defined on mesh vertices, computes the total length of the curve where the function equals the specified level value(s).

Parameters:

vfuncnp.ndarray: Scalar function values at vertices, shape (n_vertices,). Must be 1D.
levelfloat or np.ndarray: Level set value(s). Can be a single float or array of level values.

Returns:

float or np.ndarray: Length of level set. Returns float for single level, array for multiple levels.

Raises:

ValueError: If vfunc is not 1-dimensional. If no lengths could be computed.

Parameters:

vfunc (ndarray)
level (float | ndarray)

Return type:

float | ndarray

level_path(vfunc, level, get_tria_idx=False, get_edges=False, n_points=None)[source]¶

Extract a single level set path as an array of 3D points.

Thin wrapper around extract_level_paths() for the common case of a single connected level curve. Raises ValueError if the level set has more (or fewer) than one connected component; use extract_level_paths() directly for those cases.

For closed loops the returned path has a duplicate last point equal to the first, so closure can be detected via np.allclose(path[0], path[-1]).

Parameters:

vfuncnp.ndarray: Scalar function values at vertices, shape (n_vertices,). Must be 1D.
levelfloat: Level set value to extract.
get_tria_idxbool, default=False: If True, also return triangle indices for each path segment.
get_edgesbool, default=False: If True, also return mesh vertex indices and barycentric positions for each point on the path.
n_pointsint or None, default=None: If given, resample the path to this many equidistant points. Cannot be combined with get_tria_idx=True or get_edges=True.

Returns:

pathnp.ndarray, shape (n, 3): 3D coordinates along the level curve. For closed loops the last point duplicates the first.
lengthfloat: Total arc length of the path.
tria_idxnp.ndarray, shape (n-1,): Triangle index for each segment. Only returned if get_tria_idx=True.
edges_vidxsnp.ndarray, shape (n, 2): Mesh vertex indices (i, j) of the intersected edge for each point. Only returned if get_edges=True.
edges_relposnp.ndarray, shape (n,): Barycentric position along each edge (0 = vertex i, 1 = vertex j). The 3D point equals (1 - t) * v_i + t * v_j. Only returned if get_edges=True.

Raises:

ValueError: If vfunc is not 1-dimensional. If the level set has more or fewer than one connected component. If n_points is combined with get_tria_idx=True or get_edges=True.

Parameters:

vfunc (ndarray)
level (float)
get_tria_idx (bool)
get_edges (bool)
n_points (int | None)

Return type:

tuple[ndarray, …]

See also

extract_level_paths: Returns all connected components as Polygon objects.

map_tfunc_to_vfunc(tfunc, weighted=False)[source]¶

Map triangle function to vertices by attributing 1/3 to each vertex.

Distributes triangle values to their three vertices. Useful for converting per-triangle data to per-vertex data.

Parameters:

tfuncnp.ndarray: Float vector or matrix of shape (n_triangles,) or (n_triangles, N) with values at triangles.
weightedbool, default=False: If False, weigh only by 1/3 (e.g., to compute vertex areas from triangle areas). If True, weigh by triangle area / 3 (e.g., to integrate a function defined on the triangles).

Returns:

np.ndarray: Function on vertices, shape (n_vertices,) or (n_vertices, N).

Raises:

ValueError: If tfunc length doesn’t match number of triangles.

Parameters:

tfunc (ndarray)
weighted (bool)

Return type:

ndarray

map_vfunc_to_tfunc(vfunc)[source]¶

Map vertex function to triangles by averaging vertex values.

Computes triangle values by averaging the values at the three vertices of each triangle.

Parameters:

vfuncnp.ndarray: Float vector or matrix of shape (n_vertices,) or (n_vertices, N) with values at vertices.

Returns:

np.ndarray: Function on triangles, shape (n_triangles,) or (n_triangles, N).

Raises:

ValueError: If vfunc length doesn’t match number of vertices.

Parameters:

vfunc (ndarray)

Return type:

ndarray

normal_offset_(d)[source]¶

Move vertices along their normals by distance d.

Displaces each vertex along its vertex normal. Useful for creating offset surfaces or inflating/deflating meshes.

Parameters:

dfloat or np.ndarray: Move distance. Can be a scalar (same distance for all vertices) or array of shape (n_vertices,) for per-vertex distances.

Raises:

ValueError: If mesh is not oriented (required for vertex normals).

Parameters:

d (float)

Return type:

None

Notes

This modifies vertices in-place without re-initializing adjacency matrices.

normalize_()[source]¶

Normalize TriaMesh to unit surface area and centroid at the origin.

Modifies the vertices in-place by: 1. Translating centroid to origin 2. Scaling to unit surface area

Raises:

ValueError: If mesh surface area is not positive.

Return type:

None

orient_()[source]¶

Re-orient triangles of manifold mesh to be consistent.

Re-orients triangles so that vertices are listed counter-clockwise when looking from above (outside). Uses a flood-fill algorithm to propagate orientation from the first triangle.

Algorithm:

Construct list for each half-edge with its triangle and edge direction
Drop boundary half-edges and find half-edge pairs
Construct sparse matrix with triangle neighbors, with entry 1 for opposite half edges and -1 for parallel half-edges (normal flip across this edge)
Flood mesh from first triangle using neighbor matrix, tracking sign
Negative sign for a triangle indicates it needs to be flipped
If global volume is negative, flip everything (first triangle was wrong)

Returns:

int: Number of triangles flipped.

Return type:

int

Notes

This modifies the triangle array in-place to achieve consistent orientation.

refine_(it=1)[source]¶

Refine the triangle mesh by placing new vertex on each edge midpoint.

Creates 4 similar triangles from each parent triangle (1-to-4 subdivision). Modifies mesh in-place.

Parameters:

itint, default=1: Number of refinement iterations. Each iteration quadruples the number of triangles.

Parameters:

it (int)

Return type:

None

rm_free_vertices_()[source]¶

Remove unused (free) vertices.

Free vertices are vertices that are not referenced by any triangle. They can produce problems when constructing, e.g., Laplace matrices.

Modifies the mesh in-place by removing unused vertices and re-indexing triangles.

Returns:

vkeepnp.ndarray: Indices (from original list) of kept vertices, shape (n_kept,).
vdelnp.ndarray: Indices of deleted (unused) vertices, shape (n_deleted,).

Raises:

ValueError: If triangle indices exceed number of vertices.

Return type:

tuple[ndarray, ndarray]

smooth_(n=1)[source]¶

Smooth the mesh iteratively in-place using Taubin smoothing.

Parameters:

nint, default=1: Number of smoothing iterations.

Parameters:

n (int)

Return type:

None

smooth_laplace(vfunc=None, n=1, lambda_=0.5, mat=None)[source]¶

Smooth the mesh or a vertex function using Laplace smoothing.

Applies iterative smoothing: v_new = (1-lambda)*v + lambda * M*v where M is the vertex-area weighted adjacency matrix.

Parameters:

vfuncnp.ndarray or None, default=None: Float vector of values at vertices, shape (n_vertices,) or (n_vertices, N). If None, uses vertex coordinates for mesh smoothing.
nint, default=1: Number of smoothing iterations.
lambda_float, default=0.5: Diffusion speed parameter in range [0, 1]. lambda_=1 reduces to weighted average of neighboring vertices, while smaller values preserve more of the original values.
matscipy.sparse.csc_matrix or None, default=None: Precomputed smoothing matrix. If None, constructs it internally.

Returns:

np.ndarray: Smoothed vertex function, same shape as input vfunc.

Raises:

ValueError: If vfunc length doesn’t match number of vertices.

Parameters:

vfunc (ndarray | None)
n (int)
lambda_ (float)
mat (csc_matrix | None)

Return type:

ndarray

smooth_taubin(vfunc=None, n=1, lambda_=0.5, mu=-0.53)[source]¶

Smooth the mesh or a vertex function using Taubin smoothing.

Taubin smoothing alternates between shrinking (positive lambda) and expanding (negative mu) steps to reduce shrinkage while smoothing. This is a low-pass filter that better preserves mesh volume.

Parameters:

vfuncnp.ndarray or None, default=None: Float vector of values at vertices, shape (n_vertices,) or (n_vertices, N). If None, uses vertex coordinates for mesh smoothing.
nint, default=1: Number of smoothing iterations (each iteration applies both lambda and mu steps).
lambda_float, default=0.5: Positive diffusion parameter for shrinking step.
mufloat, default=-0.53: Negative diffusion parameter for expanding step. Should be negative and satisfy mu < -lambda to avoid instability.

Returns:

np.ndarray: Smoothed vertex function, same shape as input vfunc.

Raises:

ValueError: If vfunc length doesn’t match number of vertices.

Parameters:

vfunc (ndarray | None)
n (int)
lambda_ (float)
mu (float)

Return type:

ndarray

References

Gabriel Taubin, “A Signal Processing Approach to Fair Surface Design”, SIGGRAPH 1995.

smooth_vfunc(vfunc, n=1)[source]¶

Smooth the mesh or a vertex function iteratively.

This is the most simple smoothing approach of a weighted average of the neighbors.

Parameters:

vfuncarray: Float vector of values at vertices, if empty, use vertex coordinates.
nint, default=1: Number of iterations for smoothing.

Returns:

vfuncarray: Smoothed surface vertex function.

tria_areas()[source]¶

Compute the area of triangles using Heron’s formula.

Heron’s formula computes the area of a triangle by using the three edge lengths.

Returns:

np.ndarray: Array with areas of each triangle, shape (n_triangles,).

Return type:

ndarray

tria_normals()[source]¶

Compute triangle normals.

Triangle normals are computed using the cross product of two edges. Ordering of triangles is important: normals point outward for counter-clockwise oriented triangles when looking from above.

Returns:

np.ndarray: Triangle normals of shape (n_triangles, 3). Each normal is normalized to unit length.

Return type:

ndarray

tria_qualities()[source]¶

Compute triangle quality for each triangle in mesh.

Quality measure: q = 4 sqrt(3) A / (e1^2 + e2^2 + e3^2) where A is the triangle area and ei the edge length of the three edges. Constants are chosen so that q=1 for the equilateral triangle.

Note

This measure is used by FEMLAB and can also be found in: R.E. Bank, PLTMG …, Frontiers in Appl. Math. (7), 1990.

Returns:

np.ndarray: Array with triangle qualities, shape (n_triangles,). Values range from 0 (degenerate) to 1 (equilateral).

Return type:

ndarray

vertex_areas()[source]¶

Compute the area associated to each vertex (1/3 of one-ring trias).

Returns:

vareasarray: Array of vertex areas.

vertex_degrees()[source]¶

Compute the vertex degrees (number of edges at each vertex).

Returns:

np.ndarray: Array of vertex degrees, shape (n_vertices,).

Return type:

ndarray

vertex_normals()[source]¶

Compute vertex normals.

Vertex normals are computed by averaging triangle normals around each vertex, weighted by the angle that each triangle contributes at that vertex. The mesh must be oriented for meaningful results.

Returns:

np.ndarray: Vertex normals of shape (n_vertices, 3). Each normal is normalized to unit length.

Raises:

ValueError: If mesh is not oriented.

Return type:

ndarray

volume()[source]¶

Compute the volume of closed triangle mesh, summing tetrahedra at origin.

The mesh must be closed (no boundary edges) and oriented (consistent triangle orientation).

Returns:

float: Total enclosed volume.

Raises:

ValueError: If mesh is not closed or not oriented.

Return type:

float

write_fssurf(filename, image=None, coords_are_voxels=None)[source]¶

Save as Freesurfer Surface Geometry file (wrap Nibabel).

Parameters:

filenamestr: Filename to save to.
imagestr, object, None: Path to image, nibabel image object, or image header. If specified, volume_info will be extracted from the image header, and by default, vertices are assumed to be in voxel coordinates and will be converted to surface RAS (tkr) before saving. The expected order of coordinates is (x, y, z) matching the image voxel indices in nibabel.
coords_are_voxelsbool or None, default=None: Specifies whether vertices are in voxel coordinates. If None (default), the behavior is inferred: when image is provided, vertices are assumed to be in voxel space and converted to surface RAS; when image is not provided, vertices are assumed to already be in surface RAS. Set explicitly to True to force conversion (requires image), or False to skip conversion even when image is provided (only extracts volume_info).

Raises:

ValueError: If coords_are_voxels is explicitly True but image is None. If image header cannot be processed.
IOError: If file cannot be written.

Parameters:

filename (str)
image (object | None)
coords_are_voxels (bool | None)

Return type:

None

Notes

The surface RAS (tkr) transform is obtained from a header that implements get_vox2ras_tkr() (e.g., MGHHeader). For other header types (NIfTI1/2, Analyze/SPM, etc.), we attempt conversion via MGHHeader.from_header.

write_gifti(filename)[source]¶

Save mesh as a GIFTI surface file.

Writes the mesh as a .surf.gii GIFTI file using nibabel.

Parameters:

filenamestr: Output filename (conventionally ends with .surf.gii).

Raises:

ImportError: If nibabel is not installed.
OSError: If the file cannot be written.

Parameters:

filename (str)

Return type:

None

Notes

Vertex coordinates are stored as float32 and triangle indices as int32, matching the conventions of most neuroimaging software.

Examples

>>> from lapy import TriaMesh
>>> tria = TriaMesh.read_off("my_mesh.off")
>>> tria.write_gifti("my_mesh.surf.gii")

write_vtk(filename)[source]¶

Save mesh as VTK file.

Parameters:

filenamestr: Filename to save to.

Raises:

IOError: If file cannot be written.

Parameters:

filename (str)

Return type:

None