Compute Harmonic weights on a triangular mesh

Here I describe the discreet Laplace-Beltrami operator for a triangle mesh and how to derive the Laplacian matrix for that mesh. Then I provide [  C++ code ] to compute harmonic weights over a triangular mesh by solving the Laplace equation.

Motivation

The goal here is to present the Laplace operator with a concrete application, that is, computing harmonic weights over a triangle mesh. In geometry processing the Laplace operator \( \Delta \) is a extremely powerful tool. For instance it can be handy to compute or process weight maps.

A weight map is a function \( f(u, v) \) that associates a point of a surface to a real value \( f(u, v): \mathbb R^2 \rightarrow \mathbb R\). Usually the weight map is discrete, meaning if our surface is represented with triangles and vertices we can associate each vertex to a single scalar value (\(f_0=0.1\), \(f_1=-0.3\) etc ):

They are numerous occasion where you would want to associate some weight to your vertices, like setting a mass to each vertex for physic simulations, a height to create wrinkles like in bump maps, or joint influence strength like for skin weights etc. We can make use of the Laplace operator \( \Delta \) to generate weight maps \( f \). We can also smooth or process them in many ways. Here I present how to generate bounded harmonic weights by solving the Laplace equation.

Discrete formulation (per vertex)

Let's remember the discrete formulation of the Laplace operator over a triangle mesh at a single vertex \( i \):

$$ \Delta f( \vec v_i) =  \frac{1} {2A_i} \sum\limits_{{j \in \mathcal N(i)}} { (\cot \alpha_{ij} + \cot \beta_{ij}) (f_j - f_i) } $$

Barycentric cell area: Probably the simplest way to compute \( A_i \) just sum the triangles' areas adjacent to \(i\) times 3:

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

Mixed Voronoi cell area: more robust but a bit more painful to compute

float cell_area = 0 // A_i
for( each triangle tri_j from the 1-ring neighborhood of vertex i )
    if( tri_j is non-obtuse) // Voronoi safe
        cell_area += voronoi region of i in tri_j 
    else // Voronoi inappropriate
        if( the angle of 'tri_j' at 'i' is obtuse) 
            cell_area += area(tri_j)/2.0
        else 
            cell_area += area(tri_j)/4.0

 The Voronoi region of a triangle with vertex \(\vec p\),\(\vec q\) and \(\vec r\):

$$\frac{1}{8} (\|\vec p - \vec r\|^2 \cot( \text{angle at } \vec q) + \|\vec p - \vec q\|^2 \cot( \text{angle at } \vec r) )$$

Matrix representation (development)

Here things get interesting, we will represent \( \Delta f( \vec v_i) \) at every vertices using matrices so that \( \Delta f = \mathbf A \mathbf x\) with \( \mathbf x = \{ f_0, f_1 \ldots f_i \ldots \} \). We will be able to re-use this matrix in many FEM applications (Finite Element Method) and solve awesome PDEs (Partial Differential Equations) For brevity let's define:

We start from the per vertex definition refactor the expression:

$$
\begin{array}{lll}
\Delta f( \vec v_i) & = & \frac{1} {A_i} \sum\limits_{{j \in \mathcal N(i)}} { w_{ij} (f_j - f_i) }  \\
                          & = & \frac{1} {A_i} \sum\limits_{{j \in \mathcal N(i)}} { w_{ij} (f_j) } - \frac{1} {A_i}
\left ( \sum\limits_{{j \in \mathcal N(i)}} { w_{ij} } \right ) (f_i) \\
                          & = & \frac{1} {A_i} \sum\limits_{{j \in \mathcal N(i)}} { w_{ij} (f_j) } - \frac{1} {A_i}
w_{sum_i} (f_i) \\
\end{array}
$$

Now that we properly isolated our \( f_i \) and \( f_j\) (really our \(\mathbf x\)) we can infer the general matrix representation for each vertices (each matrix has the same dimension \( \mathbb R^{n \times n} \).):

$$
\mathbf A \mathbf x =
    \overbrace{ {\mathbf M}^{-1} }^{\frac{1} {A_i}} \quad
    \overbrace{ \mathbf L_{w_{ij}} \mathbf x }^{ \sum\limits_{{j \in \mathcal N(i)}} w_{ij} f_j }    
    - {\mathbf M}^{-1}
    \overbrace{ \mathbf L_{w_{sum}} \mathbf x }^{ \sum\limits_{{j \in \mathcal N(i)}} w_{ij} f_i }
$$
A example of expansion:
$$
\begin{bmatrix}
A_0 & 0        & 0 \\
0      & A_1 & 0 \\
0      & 0        & A_2\\
\end{bmatrix}^{-1}
\begin{bmatrix}
0          & w_{01} & w_{02} \\
w_{10} & 0          & w_{12} \\
w_{20} & w_{21} & 0          \\
\end{bmatrix}
\begin{bmatrix}
f_0 \\
f_1 \\
f_2 \\
\end{bmatrix}
-
\begin{bmatrix}
A_0 & 0        & 0 \\
0      & A_1 & 0 \\
0      & 0        & A_2\\
\end{bmatrix}^{-1}
\begin{bmatrix}
w_{sum_0} & 0        & 0 \\
0      & w_{sum_1} & 0 \\
0      & 0        &w_{sum_2}\\
\end{bmatrix}
\begin{bmatrix}
f_0 \\
f_1 \\
f_2 \\
\end{bmatrix}
$$

The Mass matrix \( \mathbf M \) is a diagonal matrix that stores the cell area of each vertex; elements outside the diagonal are null:

$$
\mathbf M =
\begin{bmatrix}
A_0 & 0        & 0 \\
0      & \ddots & 0 \\
0      & 0        & A_n\\
\end{bmatrix}
$$

The matrix \( \mathbf L_{w_{sum}} \) is also a diagonal matrix that stores the sum of the cotangent weights at each vertex \( i \):

$$
\mathbf L_{w_{sum}} =
\begin{bmatrix}
w_{sum_0} & 0        & 0 \\
0      & \ddots & 0 \\
0      & 0        &w_{sum_n}\\
\end{bmatrix}
$$

Lastly the Adjacency matrix \( \mathbf L_{w_{ij}} \) is symetric and each row \( i \) contains the list of adjacent vertex to weights \( w_{ij} \) to \( i \):



$$
\mathbf L_{w_{ij}} =
\begin{bmatrix}
0          & w_{01} & w_{02} & w_{03} & 0          & 0 \\
w_{10} & 0          & w_{12} & w_{13} & 0          & 0 \\
w_{20} & w_{21} & 0          & w_{23} & w_{24} & 0 \\
w_{30} & w_{31} & w_{32} & 0          & w_{34} & w_{35} \\
0          & 0          & w_{42} & w_{43} & 0          & w_{45} \\
0          & 0          & 0          & w_{53} & w_{54} & 0
\end{bmatrix}
$$

The Laplacian matrix definition (digest)

$$
\begin{array}{lll}
\mathbf A \mathbf x & = & {\mathbf M}^{-1} \mathbf L_{w_{ij}} \mathbf x - {\mathbf M}^{-1} \mathbf L_{w_{sum}} \mathbf x  \\
& = & ( {\mathbf M}^{-1} \mathbf L_{w_{ij}} - {\mathbf M}^{-1} \mathbf L_{w_{sum}}) \mathbf x  \\
& = & {\mathbf M}^{-1} (\mathbf L_{w_{ij}} - \mathbf L_{w_{sum}}) \mathbf x  \\
\end{array}
$$

What is commonly called the Laplacian matrix in the literature of geometry processing is:

$$\mathbf L = \mathbf L_{w_{ij}} - \mathbf L_{w_{sum}}$$

$$
\mathbf L_{ij} =

\left \{
\begin{matrix}
-w_{sum} & = & -\sum\limits_{{j \in \mathcal N(i)}} { w_{ij} } & \text{ when } i = j \\
w_{ij} & = & \frac{1}{2} \cot \alpha_{ij} + \cot \beta_{ij}      & \text{ if i adjacent to 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$$

Finally be aware that the general Laplacian matrix as defined in the graph theory (as opposed to geometry processing) is usually defined as \( -\mathbf L \)

Harmonic weights

To find the values \( \mathbf x = \{ \ldots f_i \ldots \} \) that describe an harmonic function (as linear as possible on x, y z coordinates. We must solve the Laplace equation:

$$
\begin{array}{lll}
\mathbf {\Delta f} & = & 0 \\
{\mathbf M}^{-1} \mathbf L \mathbf x & = & [ 0 ] \\
\mathbf L \mathbf x & = & [ 0 ] \\
\end{array}
$$

\( [ 0 ] \in \mathbb R^n \) is a column vector of zeros
The nice thing is since \( \mathbf L \) is symmetric by construction when we multiply by \( \mathbf M \) we make the linear system of equations simpler. \( \mathbf L \) is also a sparse matrix so we can take advantage of sparse solver to solve the system even quicker.

Boundary conditions: you need to set  a close region such as:

TODO image

To fix a value \( v \) of a specific vertex \( i \) simply set the row \( i \) to 0 and the diagonal element \( \mathbf L_{i,i} = 0\). then the in the column vector of zero we set \( v\) to the entry corresponding to the ith vertex.

Note: go here to find harmonic weights on a regularly spaced grid.

C++ Code

I provide [ code ] to solve the laplace equation with Eigen results are displayed with free GLUT. In solve_laplace_equation.cpp you will find code similar to the snippet below, that builds the Laplacian matrix with cotan weights set the boundary conditions by zeroing out rows of the final matrix setting its diagonal element to 1 and right hand side of the equation to the desired value.

If your mesh is represented with an half-edge data structure (each vertex knows its direct neighboors) 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;
    }
}

References