lapy.TriaMesh

class lapy.TriaMesh(v, t, fsinfo=None)

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.
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.
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_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 levelset of vfunc at a specific level as a path 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_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_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) to load mesh data from files.

area()

Compute the total surface area of triangle mesh.

Returns:

float

Total surface area (sum of all triangle areas).

avg_edge_length()

Compute the average edge length of the mesh.

Returns:

float

Average edge length.

boundary_loops()

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.

centroid()

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.

connected_components()

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.

construct_adj_dir_tidx()

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

curvature(smoothit=3)

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

curvature_tria(smoothit=3)

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

edges(with_boundary=False)

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.

euler()

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.

get_vertices(original_dim=False)

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.

has_free_vertices()

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.

is_2d()

Check if the mesh was created with 2D vertices.

Returns:

bool

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

is_closed()

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

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

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.

keep_largest_connected_component_(clean=True)

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.

level_length(vfunc, level)

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.

level_path(vfunc, level, get_tria_idx=False, get_edges=False, n_points=None)

Extract levelset of vfunc at a specific level as a path of 3D points.

For a given real-valued scalar map on the surface mesh (vfunc), this function computes the edges that intersect with a given level set (level). It then finds the point on each mesh edge where the level set intersects. The points are sorted and returned as an ordered array of 3D coordinates together with the length of the level set path.

Note

Only works for level sets that represent a single non-intersecting path with exactly one start and one endpoint (e.g., not closed)!

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 array of triangle indices for each path segment.

get_edgesbool, default=False

If True, also return arrays of vertex indices (i,j) and relative positions for each 3D point along the intersecting edge.

n_pointsint or None, default=None

If specified, resample level set into n equidistant points. Cannot be combined with get_tria_idx=True or get_edges=True.

Returns:

pathnp.ndarray

Array of shape (n, 3) containing 3D coordinates of vertices on the level path.

lengthfloat

Length of the level set.

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

edges_vidxsnp.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 level set at this point. Only returned if get_edges is True.

edges_relposnp.ndarray

Array of floats defining the relative position for each 3D point along the edges of the original mesh (defined by edges_vidxs). Value 0 corresponds to first vertex, 1 to second vertex. The 3D point is computed as: (1 - relpos) * v_i + relpos * v_j. Only returned if get_edges is True.

Raises:

ValueError

If vfunc is not 1-dimensional. If n_points is combined with get_tria_idx=True. If n_points is combined with get_edges=True.

map_tfunc_to_vfunc(tfunc, weighted=False)

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.

map_vfunc_to_tfunc(vfunc)

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.

normal_offset_(d)

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

Notes

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

normalize_()

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.

orient_()

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:

1. Construct list for each half-edge with its triangle and edge direction

2. Drop boundary half-edges and find half-edge pairs

3. Construct sparse matrix with triangle neighbors, with entry 1 for opposite half edges and -1 for parallel half-edges (normal flip across this edge)

4. Flood mesh from first triangle using neighbor matrix, tracking sign

5. Negative sign for a triangle indicates it needs to be flipped

6. If global volume is negative, flip everything (first triangle was wrong)

Returns:

int

Number of triangles flipped.

Notes

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

classmethod read_fssurf(filename)

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_off(filename)

Load triangle mesh from OFF txt file.

Parameters:

filenamestr

Filename to load.

Returns:

TriaMesh

Loaded triangle mesh.

classmethod read_vtk(filename)

Load triangle mesh from VTK txt file.

Parameters:

filenamestr

Filename to load.

Returns:

TriaMesh

Loaded triangle mesh.

refine_(it=1)

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.

rm_free_vertices_()

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.

smooth_(n=1)

Smooth the mesh iteratively in-place using Taubin smoothing.

Parameters:

nint, default=1

Number of smoothing iterations.

smooth_laplace(vfunc=None, n=1, lambda_=0.5, mat=None)

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.

smooth_taubin(vfunc=None, n=1, lambda_=0.5, mu=-0.53)

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.

References

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

smooth_vfunc(vfunc, n=1)

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

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

tria_normals()

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.

tria_qualities()

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

vertex_areas()

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

Returns:

vareasarray

Array of vertex areas.

vertex_degrees()

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

Returns:

np.ndarray

Array of vertex degrees, shape (n_vertices,).

vertex_normals()

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.

volume()

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.

write_fssurf(filename, image=None)

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, the vertices are assumed to be in voxel coordinates and are converted to surface RAS (tkr) coordinates before saving. The expected order of coordinates is (x, y, z) matching the image voxel indices in nibabel.

Raises:

IOError

If file cannot be written.

ValueError

If image header cannot be processed.

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_vtk(filename)

Save mesh as VTK file.

Parameters:

filenamestr

Filename to save to.

Raises:

IOError

If file cannot be written.