TriaMesh Geodesics

from lapy import TriaMesh

Load a triangle mesh of a flat square (OFF file).

T = TriaMesh.read_off("../data/")
# import plotting functions
import plotly
# plotly.offline.init_notebook_mode(connected=True)
from lapy.plot import plot_tria_mesh
import as pio
pio.renderers.default = "sphinx_gallery"

We now plot the triangle mesh with a function overlay of the triangle quality. Note, that this is a function attached to the triangles, not the vertices, so it is piecewise flat.

q = T.tria_qualities()
plot_tria_mesh(T,plot_edges=True, tfunc=q)


Next, we check a few properties of eigenfunctions. For this we get the first few and also the stiffness matrix (A) and the lumped mass matrix (B, reduced to a diagonal), of which we can easily compute the inverse.

# compute first eigenfunction
from lapy import Solver
fem = Solver(T,lump=True)
eval, evec = fem.eigs()
vfunc = evec[:,1]

# also get A,B (lumped), and inverse of B (easy as it is diagonal due to lumping)
A, B = fem.stiffness, fem.mass
Bi = B.copy() **= -1
TriaMesh with regular Laplace-Beltrami
Solver: spsolve (LU decomposition) ...

The mass matrix B represents the inner product so that the integral of the product of two functions x and y over the whole mesh is x B y’. The lumped mass matrix that we use here is a simplified version, as all off-diagonal elements are added to the diagonal. The entries on the diagonal represent the vertex areas and their sum is the total area of the mesh. For our unit square it will be 1.


Let’s see what happens when we apply the Laplace operator to an eigenfunction. Eigenfunctions are solutions to Delta f = - lambda f, so we should obtain a scaled version of f.


This is the same as the corresponding eigenvalue times the eigenfunction.


Laplace is also defined as the -div(grad(f)). So first applying the gradient and then the divergence to an eigenfunction and then multiplying with inv(B) should yield the same result as above again. Note, that multiplying with inv(B) is necessary to get back from the integrated divergence to the original function.

from lapy.diffgeo import compute_gradient
from lapy.diffgeo import compute_divergence
grad = compute_gradient(T,vfunc)
divx = -compute_divergence(T,grad)


Now we will replicate the idea of geodesics in heat, where first a heat diffusion is solved and massaged in the right way to yield an approximation of geodesics on the mesh. This also works on curved meshes, but for simplicity we keep using the square here. So let’s start with computing the heat diffusion from the boundary (with default time factor m=1).

from lapy import heat
bvert = T.boundary_loops()
u = heat.diffusion(T,bvert,m=1)
TriaMesh with regular Laplace-Beltrami
Matrix Format now:  csc
Solver: spsolve (LU decomposition) ...

We show some of the level sets. Note, that they are not evenly spaced and get steeper closer to the boundary.


Next step is to compute the gradient (vector field) of the heat diffusion function and normalize all vectors to unit length.

import numpy as np

# compute gradient of heat diffusion
tfunc = compute_gradient(T,u)

# normalize gradient
X = -tfunc / np.sqrt((tfunc**2).sum(1))[:,np.newaxis]
X = np.nan_to_num(X)

Then we get the integrated divergence of the normalized gradients.

divx = compute_divergence(T,X)

Finally, to obtain the distance function, we need to solve a Poisson equation. The solution can be shifted arbitrary, so we need to subtract the minimum, which should be along the boundary of the mesh.

# compute distance
from scipy.sparse.linalg import splu
useCholmod = True
    from sksparse.cholmod import cholesky
except ImportError:
    useCholmod = False

fem = Solver(T,lump=True)
A, B = fem.stiffness, fem.mass


# solve H x = b0
# we don't need the B matrix here, as divx is the integrated divergence
print("Matrix Format now: "+H.getformat())
if useCholmod:
    print("Solver: cholesky decomp - performance optimal ...")
    chol = cholesky(H)
    x = chol(b0)
    print("Solver: spsolve (LU decomp) - performance not optimal ...")
    lu = splu(H)
    x = lu.solve(b0)

# remove shift
x = x-min(x)
TriaMesh with regular Laplace-Beltrami
Matrix Format now: csc
Solver: spsolve (LU decomp) - performance not optimal ...

In short, the idea is to first get a function (heat diffusion) that flows from a point or from a set of points (boundary) through the mesh, then normalize its gradient, compute the divergence and finally step backward through the Laplace operator to find a function that has this normalized gradient. If we look at it, we actually notice that level sets are equally spaced now.