Thanks so much for your help Daniel! The code you kindly share has helped me so many projects, and putting in work to understand your detailed explanations taught me a lot.
I found the cotangent laplacian code from PyTorch3D and implemented a Numpy version below.
A big question-mark for this code is the choice of [1, 2, 0] and [2, 0, 1]. I think this is just to get the index-pairs of the edges but, if I replace the former with [0, 1, 2], the resulting laplacian will change.
The linked pytorch3d.ops.laplacian_matrices — PyTorch3D documentation PyTorch3D is a ML library, but the code itself doesn’t have much to do with ML (besides being differentiable, the documentation says)
I’m not sure how useful (fast) this is for a GH project, but learning some math along the way was fun.
from numpy.typing import NDArray
def cotangent_laplacian(vertices: NDArray, faces: NDArray, epsilon: float=1e-12):
vertices_count, faces_count = vertices.shape[0], faces.shape[0]
vertices_grouped_by_faces = vertices[faces]
v0, v1, v2 = vertices_grouped_by_faces[:, 0], vertices_grouped_by_faces[:, 1], vertices_grouped_by_faces[:, 2]
# use vector norm to find length of face edges
A : NDArray= np.linalg.norm((v1 - v2), axis=1)
B : NDArray= np.linalg.norm((v0 - v2), axis=1)
C : NDArray= np.linalg.norm((v0 - v1), axis=1)
'''Heron's formula'''
s = 0.5 * (A + B + C)
area = np.sqrt((s * (s - A) * (s - B) * (s - C)).clip(min=epsilon))
'''Law of cosines, in cotangnet'''
A_squared, B_squared, C_squared = A * A, B * B, C * C
cotangent_a = (B_squared + C_squared - A_squared) / area
cotangent_b = (A_squared + C_squared - B_squared) / area
cotangent_c = (A_squared + B_squared - C_squared) / area
cot = np.stack([cotangent_a, cotangent_b, cotangent_c], axis=1)
cot /= 4.0
ii = faces[:, [1, 2, 0]]
jj = faces[:, [2, 0, 1]]
idx = np.stack([ii, jj], axis=0).reshape((2, faces_count * 3))
m = np.zeros(shape=(vertices_count, vertices_count))
m[idx[0], idx[1]] = cot.reshape(-1)
m += m.T
return m
m = cotangent_laplacian(vertices, faces)
m.shape
