lapy.diffgeo

Differential Geometry Functions for meshes.

This module includes gradient, divergence, curvature, geodesics, mean curvature flow etc.

Note, the interface is not yet final, some functions are split into tet and tria versions.

lapy.diffgeo.compute_divergence(geom, vfunc)[source]

Compute divergence of a vertex function f.

Parameters:
geomTriaMesh or TetMesh

Mesh geometry.

vfuncnp.ndarray

3D vector field on elements, shape (n_elements, 3) or (n_elements, n_functions, 3).

Returns:
np.ndarray

Scalar function of divergence at vertices. Shape is (n_vertices,) for a single vector field or (n_vertices, n_functions) for multiple vector fields.

Raises:
ValueError

If unknown geometry type.

Parameters:
Return type:

ndarray

lapy.diffgeo.compute_geodesic_f(geom, vfunc, use_cholmod=False)[source]

Compute function with normalized gradient (geodesic distance).

Computes gradient, normalizes it, and computes function with this normalized gradient by solving the Poisson equation with the divergence of grad. This idea is also described in the paper “Geodesics in Heat” for triangles.

Parameters:
geomTriaMesh or TetMesh

Mesh geometry.

vfuncnp.ndarray

Scalar function at vertices, shape (n_vertices,) or (n_vertices, n_functions).

use_cholmodbool, default=False

If True, use the sparse Cholesky decomposition of CHOLMOD (requires cholespy) for the linear solve. If False, use the LU decomposition.

Returns:
np.ndarray

Scalar geodesic function at vertices, shape (n_vertices,) or (n_vertices, n_functions).

Parameters:
Return type:

ndarray

lapy.diffgeo.compute_gradient(geom, vfunc)[source]

Compute gradient of a vertex function f.

Parameters:
geomTriaMesh or TetMesh

Mesh geometry.

vfuncnp.ndarray

Scalar function at vertices, shape (n_vertices,) or (n_vertices, n_functions).

Returns:
np.ndarray

3D gradient vector at each element. Shape is (n_elements, 3) for a single function or (n_elements, n_functions, 3) for multiple functions.

Raises:
ValueError

If unknown geometry type.

Parameters:
Return type:

ndarray

lapy.diffgeo.compute_rotated_f(geom, vfunc)[source]

Compute function whose level sets are orthgonal to the ones of vfunc.

Parameters:
geomTriaMesh

Mesh geometry.

vfuncnp.ndarray

Scalar function at vertices, shape (n_vertices,) or (n_vertices, n_functions).

Returns:
np.ndarray

Rotated function at vertices, shape (n_vertices,) or (n_vertices, n_functions).

Raises:
ValueError

If unknown geometry type.

Parameters:
Return type:

ndarray

lapy.diffgeo.flipped_triangles(sphere, t)[source]

Find the triangles of a map onto a sphere that face its centre.

For a map of a closed surface onto a sphere centred at the origin, whose triangles keep the orientation of the surface, these are the triangles where the map folds over.

Parameters:
spherenp.ndarray

Vertex positions on a sphere centred at the origin, shape (n_vertices, 3).

tnp.ndarray

Triangle array of shape (n_triangles, 3).

Returns:
np.ndarray

Boolean mask of shape (n_triangles,), True for triangles whose normal points towards the centre.

Parameters:
Return type:

ndarray

lapy.diffgeo.spherical_spectral_map(tria, flow_iter=3, tol=1e-06, use_cholmod=False)[source]

Map a closed surface onto the unit sphere with its first three eigenfunctions.

The first three non-constant Laplace-Beltrami eigenfunctions, each scaled to [-1, 1], serve as coordinates. A few steps of mean curvature flow make this embedding more convex, and every vertex is then projected radially onto the unit sphere. The map is mirrored if needed so that most triangles keep their orientation; triangles where it folds over can be found with flipped_triangles(). Unlike tria_spherical_project(), it assumes no particular orientation of the surface in space and does not reject maps that fold.

Parameters:
triaTriaMesh

Closed triangle mesh.

flow_iterint, default=3

Mean curvature flow steps on the embedding.

tolfloat, default=1e-6

Relative accuracy of the eigenvalues; the map needs far less than machine precision.

use_cholmodbool, default=False

If True, use the sparse Cholesky decomposition of CHOLMOD (requires cholespy).

Returns:
np.ndarray

Positions on the unit sphere, shape (n_vertices, 3).

Raises:
ValueError

If the mesh is not closed.

Parameters:
Return type:

ndarray

lapy.diffgeo.tet_compute_divergence(tet, tfunc)[source]

Compute integrated divergence of a 3D tetra function f (for each vertex).

Divergence is the flux density leaving or entering a point. It can be measured by summing the dot product of the vector field with the normals to the outer faces of the 1-ring tetras around a vertex. Summing < tfunc , n_tria_oposite_v >, this is the integrated divergence, you may want to multiply with \(B^{-1}\) to get back the function in some applications.

Parameters:
tetTetMesh

Tetrahedral mesh.

tfuncnp.ndarray

3D vector field on tetrahedra, shape (n_tetrahedra, 3) or (n_tetrahedra, n_functions, 3).

Returns:
np.ndarray

Scalar function of divergence at vertices. Shape is (n_vertices,) for a single vector field or (n_vertices, n_functions) for multiple vector fields.

Parameters:
Return type:

ndarray

lapy.diffgeo.tet_compute_gradient(tet, vfunc)[source]

Compute gradient of a vertex function f (for each tetrahedron).

For a tetrahedron \((v_i, v_j, v_k, v_h)\) with volume V we have:

\[\begin{split}grad(f) &= [ (f_j - f_i) (v_i-v_k) \times (v_h-v_k) \\ & + (f_k - f_i) (v_i-v_h) \times (v_j-v_h) \\ & + (f_h - f_i) (v_k-v_i) \times (v_j-v_i) ] / (6 V).\end{split}\]
Parameters:
tetTetMesh

Tetrahedral mesh.

vfuncnp.ndarray

Scalar function at vertices, shape (n_vertices,) or (n_vertices, n_functions).

Returns:
np.ndarray

3D gradient vector at tetrahedra. Shape is (n_tetrahedra, 3) for a single function or (n_tetrahedra, n_functions, 3) for multiple functions.

Parameters:
Return type:

ndarray

Notes

Numexpr could speed up this functions if necessary.

lapy.diffgeo.tria_compute_divergence(tria, tfunc)[source]

Compute integrated divergence of a 3D triangle function f (for each vertex).

Divergence is the flux density leaving or entering a point. Note: this is the integrated divergence, you may want to multiply with \(B^{-1}\) to get back the function in some applications

Parameters:
triaTriaMesh

Triangle mesh.

tfuncnp.ndarray

3D vector field on triangles, shape (n_triangles, 3) or (n_triangles, n_functions, 3).

Returns:
np.ndarray

Scalar function of divergence at vertices. Shape is (n_vertices,) for a single vector field or (n_vertices, n_functions) for multiple vector fields.

Parameters:
Return type:

ndarray

Notes

Numexpr could speed up this functions if necessary.

lapy.diffgeo.tria_compute_divergence2(tria, tfunc)[source]

Compute integrated divergence of a 3D triangle function f (for each vertex).

Divergence is the flux density leaving or entering a point. It can be measured by summing the dot product of the vector field with the normals to the outer edges of the 1-ring triangles around a vertex. Summing < tfunc , e_ij cross n >, this is the integrated divergence, you may want to multiply with \(B^{-1}\) to get back the function in some applications.

Parameters:
triaTriaMesh

Triangle mesh.

tfuncnp.ndarray

3D vector field on triangles, shape (n_triangles, 3) or (n_triangles, n_functions, 3).

Returns:
np.ndarray

Scalar function of divergence at vertices. Shape is (n_vertices,) for a single vector field or (n_vertices, n_functions) for multiple vector fields.

Parameters:
Return type:

ndarray

Notes

Numexpr could speed-up this functions if necessary.

lapy.diffgeo.tria_compute_geodesic_f(tria, vfunc, use_cholmod=False)[source]

Compute geodesic distance on a triangle mesh.

Follows the idea of “Geodesics in Heat” by solving the Poisson equation with the divergence of the normalized gradient to obtain the geodesic distance function.

Parameters:
triaTriaMesh

Triangle mesh geometry.

vfuncnp.ndarray

Scalar function at vertices, shape (n_vertices,) or (n_vertices, n_functions).

use_cholmodbool, default=False

If True, use the sparse Cholesky decomposition of CHOLMOD (requires cholespy) for the linear solve. If False, use the LU decomposition.

Returns:
np.ndarray

Scalar geodesic function at vertices, shape (n_vertices,) or (n_vertices, n_functions).

Parameters:
Return type:

ndarray

lapy.diffgeo.tria_compute_gradient(tria, vfunc)[source]

Compute gradient of a vertex function f (for each triangle).

\[\begin{split}grad(f) &= [ (f_j - f_i) (v_i-v_k)' + (f_k - f_i) (v_j-v_i)' ] / (2 A) \\ &= [ f_i (v_k-v_j)' + f_j (v_i-v_k)' + f_k (v_j-v_i)' ] / (2 A)\end{split}\]

for triangle \((v_i,v_j,v_k)\) with area \(A\), where (.)’ is 90 degrees rotated edge, which is equal to cross(n,vec).

Parameters:
triaTriaMesh

Triangle mesh.

vfuncnp.ndarray

Scalar function at vertices, shape (n_vertices,) or (n_vertices, n_functions).

Returns:
np.ndarray

3D gradient vector at triangles. Shape is (n_triangles, 3) for a single function or (n_triangles, n_functions, 3) for multiple functions.

Parameters:
Return type:

ndarray

Notes

Numexpr could speed up this functions if necessary. Good background to read: http://dgd.service.tu-berlin.de/wordpress/vismathws10/2012/10/17/gradient-of-scalar-functions/ Mancinelli, Livesu, Puppo, Gradient Field Estimation on Triangle Meshes http://pers.ge.imati.cnr.it/livesu/papers/MLP18/MLP18.pdf Desbrun …

lapy.diffgeo.tria_compute_rotated_f(tria, vfunc, use_cholmod=False)[source]

Compute function with rotated gradient on a triangle mesh.

Solves the Poisson equation for the divergence of the rotated gradient function of vfunc. The rotation is 90 degrees around the mesh normals.

Parameters:
triaTriaMesh

Triangle mesh geometry.

vfuncnp.ndarray

Scalar function at vertices, shape (n_vertices,) or (n_vertices, n_functions).

use_cholmodbool, default=False

If True, use the sparse Cholesky decomposition of CHOLMOD (requires cholespy) for the linear solve. If False, use the LU decomposition.

Returns:
np.ndarray

Scalar function at vertices, shape (n_vertices,) or (n_vertices, n_functions).

Parameters:
Return type:

ndarray

lapy.diffgeo.tria_mean_curvature_flow(tria, max_iter=30, stop_eps=1e-13, step=1.0, use_cholmod=False)[source]

Flow a triangle mesh along the mean curvature normal.

Iteratively flows a triangle mesh along mean curvature normal (non-singular, see Kazhdan 2012). This uses the algorithm described in Kazhdan 2012 “Can mean curvature flow be made non-singular” which uses the Laplace-Beltrami operator but keeps the stiffness matrix (A) fixed and only adjusts the mass matrix (B) during the steps. It will normalize surface area of the mesh and translate the barycenter to the origin. Closed meshes will map to the unit sphere.

Parameters:
triaTriaMesh

Triangle mesh.

max_iterint, default=30

Maximal number of steps.

stop_epsfloat, default=1e-13

Stopping threshold.

stepfloat, default=1.0

Euler step size.

use_cholmodbool, default=False

Which solver to use. If True, use the sparse Cholesky decomposition of CHOLMOD (requires cholespy). If False, use the LU decomposition.

Returns:
TriaMesh

Triangle mesh after flow.

Raises:
ImportError

If use_cholmod is True but cholespy is not installed.

Parameters:
Return type:

TriaMesh

Notes

Only the system matrix of the first step is factorized. It changes little between steps, so later steps are solved with conjugate gradients preconditioned by that factorization, with a direct solve as fallback.

lapy.diffgeo.tria_spherical_project(tria, flow_iter=3, debug=False, tol=1e-06, use_cholmod=False)[source]

Compute first 3 non-constant eigenfunctions and project spectral embedding onto sphere.

Computes the first three non-constant eigenfunctions and then projects the spectral embedding onto a sphere. This works when the first functions have a single closed zero level set, splitting the mesh into two domains each. Depending on the original shape triangles could get inverted. We also flip the functions according to the axes that they are aligned with for the special case of brain surfaces in FreeSurfer coordinates.

Parameters:
triaTriaMesh

Triangle mesh.

flow_iterint, default=3

Mean curvature flow iterations (3 should be enough).

debugbool, default=False

If True, writes debug eigenvalue file.

tolfloat, default=1e-6

Relative accuracy of the eigenvalues; 0 means machine precision.

use_cholmodbool, default=False

If True, use the sparse Cholesky decomposition of CHOLMOD (requires cholespy).

Returns:
TriaMesh

Triangle mesh projected onto sphere.

Raises:
ValueError

If mesh is not closed. If direction 1 is not anterior-posterior. If global normal flip is detected. If sphere area fraction is below 0.99. If flipped area fraction is above 0.0008. If spatial volume (orthogonality) is below 0.6.

Parameters:
Return type:

TriaMesh

See also

spherical_spectral_map

The same map without the orientation for brain surfaces and checks.