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 3 float coordinates.

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

fsinfodict | None

FreeSurfer Surface Header Info.

Attributes

v	(array_like) List of lists of 3 float coordinates.
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.

Methods

area()	Compute the total surface area of triangle mesh.
avg_edge_length()	Compute the average edge length of the mesh.
boundary_loops()	Compute a tuple of boundary loops.
centroid()	Compute centroid of triangle mesh as a weighted average of triangle centers.
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 (orthogonal and in tria plane).
edges([with_boundary])	Compute vertices and adjacent triangle ids for each edge.
euler()	Compute the Euler Characteristic.
has_free_vertices()	Check if the vertex list has more vertices than what is used in tria.
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.
level_length(vfunc, level)	Compute the length of level sets.
level_path(vfunc, level[, get_tria_idx, ...])	Extract levelset of vfund at specific level as a path of 3D points.
map_tfunc_to_vfunc(tfunc[, weighted])	Map tria function to vertices by attributing 1/3 to each vertex of triangle.
map_vfunc_to_tfunc(vfunc)	Map vertex function to triangles by attributing 1/3 to each.
normal_offset_(d)	Move vertices along normal 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.
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 where.
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)	Save as Freesurfer Surface Geometry file (wrap Nibabel).
write_vtk(filename)	Save as VTK file.

Notes

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

area()

Compute the total surface area of triangle mesh.

Returns:

float

Total surface area.

avg_edge_length()

Compute the average edge length of the mesh.

Returns:

float

Average edge length.

boundary_loops()

Compute a tuple of boundary loops.

Meshes can have 0 or more boundary loops, which are cycles in the directed adjacency graph of the boundary edges. Works on trias only. Could fail if loops are connected via a single vertex (like a figure 8). That case needs debugging.

Returns:

loopslist of list

List of lists with boundary loops.

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:

centroidfloat

The centroid of the mesh.

totalareafloat

The total area of the mesh.

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.

Note

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

Parameters:

smoothitint

Smoothing iterations on vertex functions.

Returns:

u_minarray of shape (vnum, 3)

Minimal curvature directions.

u_maxarray of shape (vnum, 3)

Maximal curvature directions.

c_minarray

Minimal curvature.

c_maxarray

Maximal curvature.

c_meanarray

Mean curvature (c_min + c_max) / 2.0m.

c_gaussarray

Gauss curvature c_min * c_maxm.

normalsarray of shape (vnum, 3)

Normals.

curvature_tria(smoothit=3)

Compute min and max curvature and directions (orthogonal and in tria plane).

First we compute these values on vertices and then smooth there. Finally, they get mapped to the trias (averaging) and projected onto the triangle plane, and orthogonalized.

Parameters:

smoothitint

Number of smoothing iterations for curvature computation on vertices.

Returns:

u_minarray

Min curvature direction on triangles.

u_maxarray

Max curvature direction on triangles.

c_minarray

Min curvature on triangles.

c_maxarray

Max curvature on triangles.

edges(with_boundary=False)

Compute vertices and adjacent triangle ids for each edge.

Parameters:

with_boundarybool

Also work on boundary half edges, default ignore.

Returns:

vidsarray

Column array with starting and end vertex for each unique inner edge.

tidsarray

2 column array with triangle containing the half edge from vids[0,:] to vids [1,:] in first column and the neighboring triangle in the second column.

bdrvidsarray

If with_boundary is true: 2 column array with each boundary half-edge.

bdrtidsarray

If with_boundary is true: 1 column array with the associated triangle to each boundary edge.

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.

has_free_vertices()

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

Returns:

bool

Whether vertex list has more vertices or not.

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.

level_length(vfunc, level)

Compute the length of level sets.

Parameters:

vfuncarray

Float vector of values at vertices (here only scalar function 1D).

levelfloat | array

Level set value or array of level values.

Returns:

lengthfloat | array

Length of level set (or array of lengths).

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

Extract levelset of vfund at 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!

Additional options: get_tria_idx and get_edges when True will also return an array of triangle ids for each path segment, defining the triangle i, where the path from i to i+1 passes through (so for a path with n points, these will be n-1 triangle ids). Furthermore, an array of edge vertex indices of shape (n,2) can be obtained defining the two vertices of the intersecting edge in the original mesh for each 3D point on the path. A second array is returned, defining the relative position of the intersecting point along this edge as a float from 0 (start vertex) to 1 (end vertex). This information can, for example, be useful for interpolating a second surface map at the new path point coordinates. Neither of this information is available when n_points is used to resample the path into n_points equidistant new points as the association to edges or triangles in the original mesh is lost.

Parameters:

vfuncarray

Float vector of values at vertices (here only scalar function 1D).

levelfloat

Level set value.

get_tria_idxbool, default: False

Also return a list of triangle indices for each edge, default False.

get_edgesbool, default: False

Also return a list of two vertex indices (i,j) for each 3D point and a list of the relative position defining the 3D point along that edge (i,j) from the original mesh, default False.

n_pointsint

Resample level set into n equidistant points. Cannot be combined with get_tria_idx=True nor with get_edges=True.

Returns:

path: array

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

lengthfloat

Length of level set.

tria_idxarray

Array of triangle index for each segment on the path (length n-1 if path is length n). Will only be returned if get_tria_idx is True.

edges_vidxsarray

Array of shape (n,2) of vertex indices (i,j) for each 3D point, defining the vertices of the original mesh of the edge intersecting the level set at this point. Will only be returned if get_edges is True.

edges_relpos: array

Array of floats defining the relative position for each 3D point along the edges of the original mesh (defined by the two points in edges_vidxs). Float value 0 defines first point, and 1 defines end point. So the 3D point of the path is computed (1 - relpos) v_i + relpos v_j. Will only be returned if get_edges is True.

map_tfunc_to_vfunc(tfunc, weighted=False)

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

Uses vertices and trias.

Parameters:

tfuncarray

Float vector or matrix (#t x N) of values at vertices.

weightedbool, default=False

False, weigh only by 1/3, e.g. to compute vertex areas from tria areas True, weigh by triangle area / 3, e.g. to integrate a function defined on the trias, for example integrating the “one” function will also yield the vertex areas.

Returns:

vfuncarray

Function on vertices vector or matrix (#v x N).

map_vfunc_to_tfunc(vfunc)

Map vertex function to triangles by attributing 1/3 to each.

Uses number of vertices and trias

Parameters:

vfuncarray

Float vector or matrix (#t x N) of values at vertices.

Returns:

tfuncarray

Function on trias vector or matrix (#t x N).

normal_offset_(d)

Move vertices along normal by distance d.

normal_offset(d) moves vertices along normal by distance d

Parameters:

dint | array

Move distance.

normalize_()

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

Modifies the vertices.

orient_()

Re-orient triangles of manifold mesh to be consistent.

Re-orients triangles of manifold mesh to be consistent, so that vertices are listed counter-clockwise, when looking from above (outside).

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 tria using triangle neighbor matrix and keeping track of sign

• When flooded, negative sign for a triangle indicates it needs to be flipped

• If global volume is negative, flip everything (first tria was wrong)

Returns:

flippedint

Number of trias flipped.

classmethod read_fssurf(filename)

Load triangle mesh from FreeSurfer surface geometry file.

Parameters:

filenamestr

Filename to load.

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.

Thus creates 4 similar triangles from one parent triangle. Modifies mesh in place.

Parameters:

itint

Number of iterations.

rm_free_vertices_()

Remove unused (free) vertices.

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

Will update v and t in mesh.

Returns:

vkeeparray

Indices (from original list) of kept vertices.

vdelarray

Indices of deleted (unused) vertices.

smooth_(n=1)

Smooth the mesh iteratively in-place.

Parameters:

nint

Smoothing iterations.

smooth_vfunc(vfunc, n=1)

Smooth the mesh or a vertex function iteratively.

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:

areasarray

Array with areas of each triangle.

tria_normals()

Compute triangle normals.

Ordering of triangles is important: counterclockwise when looking.

Returns:

narray of shape (n_triangles, 3)

Triangle normals.

tria_qualities()

Compute triangle quality for each triangle in mesh where.

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:

array

Array with triangle qualities.

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:

vdegarray

Array of vertex degrees.

vertex_normals()

Compute vertex normals.

get_vertex_normals(v,t) computes vertex normals

Triangle normals around each vertex are averaged, weighted by the angle that they contribute. Ordering is important: counterclockwise when looking at the triangle from above.

Returns:

narray of shape (n_triangles, 3)

Vertex normals.

volume()

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

Returns:

volfloat

Total enclosed volume.

write_fssurf(filename)

Save as Freesurfer Surface Geometry file (wrap Nibabel).

Parameters:

filenamestr

Filename to save to.

write_vtk(filename)

Save as VTK file.

Parameters:

filenamestr

Filename to save to.