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 Cholesky decomposition from scikit-sparse cholmod for the linear solve. If False, use spsolve (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.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 Cholesky decomposition from scikit-sparse cholmod for the linear solve. If False, use spsolve (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 Cholesky decomposition from scikit-sparse cholmod for the linear solve. If False, use spsolve (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 Cholesky decomposition from scikit-sparse cholmod. If False, use spsolve (LU decomposition).

Returns:
TriaMesh

Triangle mesh after flow.

Raises:
ImportError

If use_cholmod is True but scikit-sparse is not installed.

Parameters:
Return type:

TriaMesh

Notes

Numexpr could speed up this functions if necessary.

lapy.diffgeo.tria_spherical_project(tria, flow_iter=3, debug=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.

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

Notes

Numexpr could speed up this functions if necessary.