Surface fairing, implementation detail

Question

I'm reading through Polygon Mesh Processing, section 4.3 (Surface fairing).

The reference I mentioned defined defines "Surface Fairing" as

The goal of surface fairing is to compute shapes that are as smooth as possible... in general fair surface should follow the principle of simplest shapes, the surfaces should be free by any unnecessary details or oscillation".

It continues...

This can be modeled by a suitable energy that penalizes unaesthetic behavior of the surface. A minimization of this fairness energy - subject to user defined constraints - eventually yields to the desired shape.

The energy defined is function of a surface $x : \Omega \to \mathcal{S}$ and the form is

$$ E_M(x) = \int_{\Omega} \sqrt{det(\bf{I})}dudv $$

where $\bf{I}$ is the first fundamental form, however instead of this the dirichlet energy

$$ \tilde{E}_M(x) = \int_{\Omega} \lVert x_u \rVert^2 + \lVert x_v \rVert^2 dudv $$

Is the one usually used. Higher order energies are also used, depending on how smooth is the shape we're looking for. Using the Euler-Lagrange equation (in a continuous setting) we end up with the simple PDE

$$ \Delta x = 0, $$

subject to suitable boundary conditions. Discretazion to meshes is done by means of the Laplace-beltrami operator which I'll define below.

In general it is my understanding that everything boils down to solve a problem of the form

$$ \Delta^k x = 0 $$

where $x = x(u,v)$ is our surface. If we descretize to a mesh the equation above might be rewritten as

$$ L^k x = 0 $$

Where $L^k$ is the $k^{th}$ power of the Laplacian matrix, and $L$ itself is defined using the Discrete Laplace Beltrami operator, which is something like

$$ \Delta x = \frac{1}{2A_T} \sum_{j} (cotg (\alpha_{i,j}) + cotg (\beta_{i,j}))(x_j - x_i) $$

However despite the fact it might be very simple to understand in principle at least I'm trying to figure how this would be implemented.

This is because to me the angles depends from the mesh vertices and therefore the system isn't actually linear. The other approach I was thinking is to use a gradient flow (diffusion) approach of the form

$$ \frac{\partial x}{\partial t} = -\Delta^k x $$

the above, in principle can be discretized as follows

$$ x(t + h) = x(t) - h L^k x(t) $$

And $L$ might be computed on the current known state of the surface $x(t)$, however the above is an explicit integrator which might lead to instability, and solving the issue would be done by solving

$$ (id + h L^k) x(t + h) = x(t) $$

But in such case there's the problem of $L$ which would depend from the unknown state $x(t+h)$, so the problem I'm having here is actually similar to the first one I exposed...

So the actual question is:

If the mesh $x$ is unknown so is the discretization of $\Delta^k$, how is this operator constructed then?

Update

I believe I've found the original paper, and I quote the bit I find interesting for my question:

From section 5.4.:

We can however proceed in exactly the same way, as the changes induced in a time step will be small. We simply compute the non-zero coefficients of the matrix: $I - \lambda dt K$ where $K$ represents the matrix of the curvature normals...

So I guess the idea is actually what has been suggested already in one of the comments, we compute the matrix $I - \lambda dt K$ for the current iteration (say $i$) and we keep the very same matrix for a small number of iterations $n$ after such number of iterations we update again the matrix.

However in the paper the autors confirm that the Laplace-Beltrami operator isn't linear and therefore the system isn't in general linear.

Could anyone confirm if my understanding is correct?

Answer 1

Yes, your understanding it correct. The Laplace-Beltrami depends on the current state $x$ and you have to recompute $L$ if $x$ changes. Therefore you cannot write matix of $\Delta$ without knowing the state $x$.

---

To expand on solving $\frac{\partial x}{\partial t} = -\Delta x$:

Discretization in space turns $\Delta$ to $L$, discretization in time turns $\frac{\partial x}{\partial t}$ to $\frac{x_{n+1} - x_n}{h}$ and decide if we evaluate the right hand side at time $t_n$ or $t_{n+1}$

Option 1 = Forward Euler: $$ x_{n+1} - x_n = h L_n x_n $$ where $L_n = L(x_n)$, $L$ is considered as a function of the state $x$. As you have mentioned this method can be unstable if we take too big step $h$. Therefore, we often resort to the second option.

Option 2 = Backward Euler: $$ x_{n+1} - x_n = h L_{n+1} x_{n+1} $$ Because $L_{n+1}$ depends on $x_{n+1}$ the above equation is non-linear system, we cannot solve it with a simple linear solve.

The common approach is to Taylor expand $L$ around the point $x_n$, general Taylor expansion around point $x$ is $$ L(x + \Delta x) = L(x) + \nabla L(x) \cdot \Delta x + \dots $$ therefore for $x= x_n$ and $\Delta x = x_{n+1} - x_n$ we get $$ L_{n+1} = L(x_{n+1}) = L(x_n + \Delta x) = L(x_n) + \nabla L(x_n) \cdot \Delta x + \dots $$

The simplest thing to do is to take just the first term of the Taylor expansion, i.e. we replace $L_{n+1}$ with $L_n$. This way, we end up with a linear problem $$ x_{n+1} - x_n = h L_{n} x_{n+1} $$

You can take higher order terms or solve the Backward Euler equation with non-linear solver. However, in your case that is probably unnecessary and the headache of computing the gradient of $L$ is not worth it(I might be wrong here, I do not have hands on experience with solving this equation).