Discretization of the Laplace Beltrami operator with the classical contangent weights

In this notebook, we implement the classical formula for the laplace beltrami operator. Following this formula, given a field of values $\left\{v_i\right\}_{i=0}^{n-1}$ attached to the mesh vertices $\left\{\mathbf{p_i}\right\}_{i=0}^{n-1}$, the Laplacian at $\mathbf{p_i}$ can be computed as

$$ \Delta{}u_i = \frac{1}{2}\sum_{\mathbf{p_j} \mbox{ voisin de } \mathbf{p_i}}(\cot{\alpha} + \cot{\beta})(u_j - u_i) $$

where $\alpha$ and $\beta$ are the angles facing edge $[\mathbf{p_i}, \mathbf{p_j}]$ in the two adjacent triangles.

%matplotlib inline

import math
import numpy as np
import scipy
from scipy import sparse
from scipy.sparse.linalg import eigsh, inv
import matplotlib
import matplotlib.pyplot as plt
from pygel3d import hmesh
from pygel3d.hmesh import Manifold
from pygel3d import jupyter_display as jd

jd.set_export_mode(True)

suzanne = hmesh.load("suzanne_triangles.obj")

jd.display(suzanne)

The Laplace Beltrami operator is the matrix which yields the vector of $\{\Delta{}u_i\}$ when multiplied by the vector of $\{u_i\}$. For triangle $(\mathbf{p_i}, \mathbf{p_j}, \mathbf{p_k})$ of the mesh and each edge $[\mathbf{p_i}, \mathbf{p_j}]$ of the triangle, the cotangent at angle $\mathbf{p_k}$ is computed, and its contribution added to the matrix.

def laplace_beltrami_operator(mesh):

    vsize = len(mesh.vertices())
    fsize = len(mesh.faces())

    C = sparse.lil_matrix((vsize, vsize))
    masses = np.zeros(vsize)

    for f in suzanne.faces():
        verts = list(mesh.circulate_face(f,'v'))
        area = mesh.area(f)
    
        for i,v in enumerate(verts):
            vi = verts[(i+1)%3]
            vj = verts[(i+2)%3]
            ei = mesh.positions()[vi] - mesh.positions()[v]
            ej = mesh.positions()[vj] - mesh.positions()[v]
        
            cot = ei.dot(ej) / (4*area)
            C[vi,vj] += cot
            C[vj,vi] += cot
            C[vi,vi] -= cot
            C[vj,vj] -= cot
        
            masses[v] += area
        

    M = sparse.dia_matrix((masses / 3, [0]), (vsize, vsize))
    L = sparse.csr_matrix(sparse.dia_matrix((3 / masses, [0]), (vsize,vsize)) * C)
    return L,M

Eigenvalues of the Laplace Beltrami operator

Eigenvalues of the Laplace Beltrami operator define the vibration modes of the object. These can be used to identify arms and legs in a mesh, of to compute features for shape matching. Matching using such features allows the matching of deformed objects, as long as the deformation is rigid.

L,M = laplace_beltrami_operator(suzanne)

D = sparse.csr_matrix(sparse.dia_matrix((1/M.diagonal(), [0]), L.shape) * L)

l, ev = eigsh(D, 100, which='SM')

jd.display(suzanne,data=ev[:,5], wireframe=False)

The Laplacian can also be used in simulating heat diffusion processes. In such a process, the variation of the temperature at each time step is provided by the Laplacian. Below, a single vertex of the mesh is heated, and we then use one step of backward Euler is used to simulate the temperature after some time. This can be used to compute eodesic distances for instance.

U0 = np.zeros(M.shape[0])
U0[0] = 1

t = 100
U = scipy.sparse.linalg.spsolve(M - t * L, M*U0)

jd.display(suzanne, data = U, wireframe = False)