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)

Compute divergence of a vertex function f.

Parameters:

geomTriaMesh or TetMesh

Mesh geometry.

vfuncnp.ndarray

Scalar function at vertices, shape (n_vertices,).

Returns:

np.ndarray

Scalar function of divergence at vertices, shape (n_vertices,).

Raises:

ValueError

If unknown geometry type.

lapy.diffgeo.compute_geodesic_f(geom, vfunc)

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

Returns:

np.ndarray

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

lapy.diffgeo.compute_gradient(geom, vfunc)

Compute gradient of a vertex function f.

Parameters:

geomTriaMesh or TetMesh

Mesh geometry.

vfuncnp.ndarray

Scalar function at vertices, shape (n_vertices,).

Returns:

np.ndarray

3D gradient vector at each element, shape (n_elements, 3).

Raises:

ValueError

If unknown geometry type.

lapy.diffgeo.compute_rotated_f(geom, vfunc)

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

Returns:

np.ndarray

Rotated function at vertices, shape (n_vertices,).

Raises:

ValueError

If unknown geometry type.

lapy.diffgeo.tet_compute_divergence(tet, tfunc)

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

Returns:

np.ndarray

Scalar function of divergence at vertices, shape (n_vertices,).

lapy.diffgeo.tet_compute_gradient(tet, vfunc)

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

Returns:

np.ndarray

3D gradient vector at tetrahedra, shape (n_tetrahedra, 3).

Notes

Numexpr could speed up this functions if necessary.

lapy.diffgeo.tria_compute_divergence(tria, tfunc)

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

Returns:

np.ndarray

Scalar function of divergence at vertices, shape (n_vertices,).

Notes

Numexpr could speed up this functions if necessary.

lapy.diffgeo.tria_compute_divergence2(tria, tfunc)

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

Returns:

np.ndarray

Scalar function of divergence at vertices, shape (n_vertices,).

Notes

Numexpr could speed-up this functions if necessary.

lapy.diffgeo.tria_compute_geodesic_f(tria, vfunc)

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

Parameters:

triaTriaMesh

Triangle mesh.

vfuncnp.ndarray

Scalar function at vertices, shape (n_vertices,).

Returns:

np.ndarray

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

lapy.diffgeo.tria_compute_gradient(tria, vfunc)

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

Returns:

np.ndarray

3D gradient vector at triangles, shape (n_triangles, 3).

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)

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

This is done by rotating the gradient around the normal by 90 degrees, then solving the Poisson equations with the divergence of rotated grad.

Parameters:

triaTriaMesh

Triangle mesh.

vfuncnp.ndarray

Scalar function at vertices, shape (n_vertices,).

Returns:

np.ndarray

Rotated scalar function at vertices, shape (n_vertices,).

Notes

Numexpr could speed up this functions if necessary.

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

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.

Notes

Numexpr could speed up this functions if necessary.

lapy.diffgeo.tria_spherical_project(tria, flow_iter=3, debug=False)

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.

Notes

Numexpr could speed up this functions if necessary.