In the book Polygon Mesh Processing, there is a section dedicated to fairing. In particular page 59 has this lovely figure:

enter image description here

I am specially interested in the rightmost image. According to the text this is done by fin ding the set of values of x such that of $\Delta ^ 3 x = 0$

The text claims that all that is needed in these cases is to solve

$$L^3x = 0$$

Where $L$ is a Laplacian matrix.

What I am not following is how to construct this Laplacian matrix, in particular, the text mentions that it is simpler to fix 2 rings of vertices for the case where you want the middle figure, so I assume you want to fix 3 rings for the rightmost figure.

I was hoping someone could explain to me how you actually feed in the coefficients of the fixed vertices into the system.

  • $\begingroup$ $L^2$ is just matrix multiplication of $L$ with itself: $L^2=L\cdot L$, similarly for $L^3$. In practice you can also just apply the matrix-vector products one at a time: $L^3v = L(L(Lv))$. Note that the above formulation is actually slightly wrong - the mistake dates back to Kobbelt's paper, there should be a mass matrix involved in there, i.e. $LM^{-1}LM^{-1}L$ instead of $L^3$. I will write a more detailed answer once I have the time. $\endgroup$
    – lightxbulb
    Sep 12, 2023 at 9:24
  • $\begingroup$ Btw it's a 2-ring neighbourhood, because iterating two times of a connectivity matrix of a 1-ring neghbourhood yields a 2-ring neighbourhood matrix. $\endgroup$
    – lightxbulb
    Sep 12, 2023 at 9:41
  • $\begingroup$ @lightxbulb thank you for the answer. I am still confused. I don't fully understand how to fill in L. In particular, in this case the Laplacian is being applied to the coordinates, so they are interdependent. So I don't fully understand how to linearized it. $\endgroup$
    – Makogan
    Sep 12, 2023 at 19:12
  • $\begingroup$ You can use the cotan formula on a mesh to compute $L$: en.wikipedia.org/wiki/Discrete_Laplace_operator#Mesh_Laplacians The mass matrix is not in general the identity, you can build a lumped mass matrix by summing the triangle areas around a vertex. As far as $x$ being position and $L(x)$ also being a function of position - this is correct. One way to solve it is to freeze $L(x^k)$ at iteration $k$ and solve for $x^{k+1}$: $L(x^k)x^{k+1} = 0$. I don't know whether something additional needs to be done here for stability. Note that the Dirichlet conditions matter. $\endgroup$
    – lightxbulb
    Sep 12, 2023 at 21:19
  • $\begingroup$ @lightxbulb sorry but this makes no sense to me. The blue area in the picture is going to vary widely depending on what one wants the final shape to be. But if you use the standard mesh laplacian, the coefficients of your matrix are set by the original guess of the mesh geometry. That sounds wrong, they should be updated over time, no? $\endgroup$
    – Makogan
    Sep 12, 2023 at 21:43

1 Answer 1


Disclaimer: I haven't read most of the papers that the book refers to, so what you will find below is my own understanding based on very little data (i.e. I may have severely misunderstood some part of the book).

Laplace-Beltrami Operator

In a finite element formulation you typically have a stiffness matrix $W$ and a mass matrix $M$. For a 2-manifold one possible discretisation is the cotangent one: $$W_{ij} = \begin{cases} -\frac{1}{2}(\cot\alpha_{ij}+\cot\beta_{ij}) &\text{for } j\in\mathcal{N}(i),\\ \sum_{k\in\mathcal{N}(i)}C_{ik} &\text{for } j=i, \\ 0 &\text{for } j\ne i \land j\notin \mathcal{N}(i), \end{cases} $$ where $\alpha_{ij}$ and $\beta_{ij}$ are the angles opposite to the edge $ij$. You can also find the same up to a factor $\frac{1}{2}$ in the appendix of the book you linked (note that in said book they use instead of $W$ and $M$, $M$ and $D$), or a few pages prior to the one you cite. In the appendix of the book that you have the mass matrix $M$ is a lumped (diagonal) matrix, where $M_{ii}$ is the Voronoi area for vertex $i$. Then finally they use $\Delta_S\approx L = -M^{-1}W$ (here $\Delta_S$ describes the Laplace-Beltrami operator on the manifold $S$).

Continuous Harmonic Formulation

As far as I can tell the continuous formulation of the fairing approach with a membrane surface would be something like this: \begin{align} \Delta_S x &= 0, \quad x\in S, \\ x &= p, \quad x\in \partial S, \end{align} where $S$ is the blue part in the images, $\partial\Omega$ is the boundary (the interface between white and blue), and $p$ are the (fixed) positions at said boundary. You may note that the above is actually a non-linear equation. This is the case because $\Delta_{S}$ depends on the surface itself, and $x$ also depends on the surface itself. Thus I would actually suggest (for the sake of devising a numerical algorithm) to instead consider the time evolution equation: \begin{alignat}{3} \partial_t x(t) &= \Delta_{S(t)} x(t), &\quad& x\in S(t), &\quad& t\in(0,\infty) \\ x(t) &= p, &\quad& x\in \partial S(t), &\quad& t\in[0,\infty), \\ S(0) &= S_0. \end{alignat} This should be equivalent to the first formulation at the steady state $\partial_t x = 0$.

Discrete Harmonic Formulation

The first continuous formulation can be discretised just as: $$(C-(I-C)L)X = CP,$$ where $\Delta_{S} \approx L = -M^{-1}W\in\mathbb{R}^{n\times n}$, $X\in\mathbb{R}^{n\times 3}$ are the $n$ vertices' positions that are inside the blue region and also the ones from the boundary (vertices that are from the white region but have a blue neighbour), $P\in\mathbb{R}^{n\times 3}$ are the initial positions. Finally $C = \operatorname{diag}(c)$, where $c \approx 1_{\partial S}$: $$c_i=\begin{cases} 1, \quad \text{$i$ is a boundary vertex}, \\ 0, \quad \text{$i$ is not a boundary vertex}. \end{cases}$$ $(C-(I-C)L)$ essentially negated $L$ and replaces its rows corresponding to Dirichlet vertices with $e_k$ (where $k$ is the index of the row being replaced). You can fully eliminate the Dirichlet vertices by considering the reduced formulation as in FEM also. Note however that in either case $L$ is actually a function of $X$ (since these cotan angles and areas depend on the position of the vertices), this means that the system is non-linear. A standard approach is to try and use a fixed point iteration, i.e. have an initial guess $X^0$, then compute $L^0$ from it, and solve the system $(C-(I-C)L^0)X^1 = CP$. Repeating this yields the iteration $(C-(I-C)L^k)X^{k+1} = CP$, and you stop when e.g. $\|X^{k+1}-X^k\|<\epsilon$ in some iteration (note that the superscripts here are iteration indices and not powers). I have no idea how stable this is however. You can of course try other approaches for solving nonlinear systems of equations.

For the time-evolution equation we can discretise the time derivative, e.g. using explicit or semi-implicit Euler. In the explicit variant you have: $$\partial_t x(t) = \Delta_{S(t)} x(t) \implies \frac{x^{k+1}-x^k}{\tau} = \Delta_{S^k} x^k,$$ while in the semi-implicit variant you have: $$\partial_t x(t) = \Delta_{S(t)} x(t) \implies \frac{x^{k+1}-x^k}{\tau} = \Delta_{S^k} x^{k+1}.$$ You can of course also try a fully-implicit variant, but that will result in a nonlinear system of equations that you must solve. Once you discretise also the spatial domain and incorporate the boundary conditions you get: \begin{alignat}{2} x^{k+1} &= (I-C)(I+\tau L^k)x^k + CP &\quad&& \text{(explicit Euler)}\\ (C+(I-C)(I-\tau L^k))x^{k+1} &= CP &\quad&& \text{(semi-implicit Euler)}. \end{alignat} The semi-implicit variant should be stable for all $\tau>0$ I believe, but larger step sizes may give poor results.

Polyharmonic Formulation

For the polyharmonic formulation you need additional boundary conditions. What I can suggest are Navier conditions: \begin{align} (-\Delta_S)^q x &= 0, \quad x\in S, \\ (-\Delta)^l x &= p^l, \quad x\in \partial S, \, 0\leq l \leq q-1. \end{align} Then this gets discretised as $$\left(C\sum_{l=0}^{q-1}(-L)^{l}+(I-C)(-L)^q\right)X = C\sum_{l=0}^{q-1}P^l.$$ Note that the superscripts here refer to powers, except for the superscript of $P$ which refers to just different matrix indices.

Implementation Details

When you solve the linear systems you do not have to solve simultaneously for each coordinate, you can solve them one at a time. The system matrix remains the same, only the right-hand side changes. You may also want symmetrize the $L$ as described in the book, then you could use a conjugate gradient solver since the matrices should be symmetric positive definite (provided you have at least 1 Dirichlet vertex).

When you have to deal with the $(-L)^j$ powers you can just implement them initially as iterated matrix-vector products. Later on you can compute the explicit form of the matrices in order to reuse cache more.

The choice of the Navier conditions for the polyharmonic case lets you directly use $(-L)^l$ (as seen in the last discretisation) instead of having to derive special modification for the boundary terms.

Also, please post some screenshots if you make this work. It's nice to know if what I mentioned actually works out and see the result of it.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.