Definition: Laplacian Matrix for triangle meshes


Quick definition of the Laplacian matrix on a triangular mesh. I give an in depth explanation here.

What is commonly called the Laplacian matrix \( \mathbf L \) in the literature of geometry processing is:

$$
\mathbf L_{ij} =
\left \{
\begin{matrix}
w_{ij} = \frac{1}{2} \cot \alpha_{ij} + \cot \beta_{ij}      & \text{ if j adjacent to i} \\
-\sum\limits_{{j \in \mathcal N(i)}} { w_{ij} } & \text{ when } i = j \\
0    & \text{ otherwise } \\
\end{matrix}
\right .
$$

On the other hand the Laplacian *operator* is defined as \( \mathbf {\Delta f} = {\mathbf M}^{-1} \mathbf L \) with the Mass matrix \( \mathbf M \), a diagonal matrix that stores the cell area (blue area on the figure) of each vertex:

$$
\mathbf M^{-1} =
\begin{bmatrix}
\frac{1}{A_0} & 0        & 0 \\
0      & \ddots & 0 \\
0      & 0        & \frac{1}{A_n}\\
\end{bmatrix}
$$

$$A_i = 3 \sum\limits_{{T_j \in \mathcal N(i)}} {area(T_j)}$$

Code

See the [ C++ code ] to build the laplacian matrix with cotan weights (get_laplacian() procedure)

If your mesh is represented with an half-edge data structure (each vertex knows its direct neighbours) the pseudo code to compute \( \mathbf L \) is:

 // angle(i,j,t) -> angle at the vertex opposite to the edge (i,j)
for(int i : vertices) {
    for(int j : one_ring(i)) {
        sum = 0;
        for(int t : triangle_on_edge(i,j))
        {
            w = cot(angle(i,j,t));
            L(i,j) = w;
            sum += w;
        }
        L(i,i) = -sum;
    }
}

On the other hand the Laplacian \( \mathbf L \) may be built by summing together contributions for each triangle, this way only the list of triangles is needed:

for(triangle t : triangles)
{
    for(edge i,j : t)
    {
        w = cot(angle(i,j,t));
        L(i,j) += w;
        L(j,i) += w;
        L(i,i) -= w;
        L(j,j) -= w;
    }
}

No comments

All html tags except <b> and <i> will be removed from your comment. You can make links by just typing the url or mail-address.
Spam bot question: