In its simplest form, it amounts to replacing a point by the average of its neighbors. For a mesh, a neighborhood can be a flexible concept, but we can reduce it to being the one-ring of each vertex, i.e., the vertices connected through an edge to the current vertex. This is valid for both triangle and quad meshes or any polygonal faces. For quad meshes, we can probably extend it to contain vertices that lie opposite (along the diagonal of the quad) to the current vertex. This is often done for cloth meshes to define stretch springs - in which case we don't really have a mesh anymore, but rather a graph. Note that this extension to quads (and cloth) is not really a second order neighborhood (vertices connected through 2 edges), but rather a subset of it.
As you can see, already the notion of a neighborhood is rich enough for a discussion. I guess we can also apply the same concept in a mesh-free context (like in SPH), where a neighborhood is defined by a radius around the current vertex. The next tricky part are the weights of each vertex in the neighborhood. In my experiments, using the "umbrella" Laplacian operator (e.g. [-1 2 1]) seems to be just fine for many cases and I'll stick with it for simplicity.
In this blog post, I will try to look at Laplacian smoothing and its various facets and try to learn something new out of it. Mainly, how is it many things at the same time? A low pass filter, a membrane simulation, a 2D diffusion process, etc. And how are they connected? Definitely, it's not going to be an exhaustive answer.
Let's start with the math:
$p'_i = p_i + \lambda L p_i = (I + \lambda L) p_i = p_i + \lambda \Delta p_i$
The new position $p'$ is the old position plus a delta term which is obtained by applying the Laplacian $L$ operator (which is in fact a matrix) to the current position $p$. But this operator depends also on the neigboring vertices of $p$. So, in the end it's a sum of the form
$\Delta p = \sum_j w_j (p_j - p_i)$
where the $w$ terms are called the weights (like in a weighted average or a convolution).
I won't talk here much about the shrinking effect of Laplacian smoothing. The most common solutions are to use the Taubin method (also known as the "lambda | mu" method) or trying to preserve the volume somehow.
Image filters
Laplacian filters are used in image processing to detect edges or sharpen the image. Usually, the convolution kernels are not "normalized" and look something like [-1 2 -1]. My guess is this is because they want to exaggerate its effect, as the image colors are clamped anyway between 0 and 255. Another important distinction is that the signed is reversed from the one used for smoothing curves and meshes. From my experiments on images and meshes so far, my conclusion is that a negative L will sharpen the field values and a positive one will smooth them.
Another interesting thing I noticed is that the resulting averaging kernel of the form (e.g. [0.5 0 0.5]) obtained after adding the Laplacian (with a full step) does not include the original value anymore. The result is still a smoothing, hence a blurring of the image, but very different from a classic Gaussian low pass filter (or sinc, or cubic B-spline for that matter). However, Taubin argues that the averaging is still a Gaussian weighting, but maybe the difference is that one works per nodes (for images) and the other per edges (for meshes). Still, it's pretty clear they are both low pass filters, so it's probably just different ways of doing it, with different derivations.
Clearly, the biggest distinction is that images are represented on a regular grid. If we were to represent an image as a scalar field broken down into triangular elements the situation would be more similar. It would be like working on the uv space of the mesh. But an image also does not have intrinsic curvature, but that is taken into account by the Laplace-Beltrami operator (the contangent weights). But if we use uniform weights, the analogy is pretty clear, even more if we used quad meshes, as then we would have something very close to a grid, at least topologically. By the way, there are two flavours of 3x3 Laplacian convolution kernels in image processing: one uses only direct neighbors and the other also the diagonal ones (as discussed above).
Mesh smoothing
The Laplacian is arising from minimizing some kind of membrane (or thin plate) energy. This in turn is usually the squared norm of some strain tensor, which also approximates total curvature. Given this objective, it is clear to see that the end result is a surface without curvature or at least the minimal curvature everywhere that still preserves the topology of the mesh. Given no other constraints, the minimization may also shrink the mesh in the process.
Desbrun et al. show that smoothing is in fact a diffusion PDE (time-dependent) integrated on a curved mesh surface. (It's like a temperature becoming uniform everywhere on the mesh; or some spiky hairs on it, like the normal vectors.) But the mesh is already a space discretization, so in the end we only need to integrate an ODE in time. This can be done either using explicit Euler resulting in the classic repeatead application of the Laplacian operator. The method using implicit Euler is called "implicit fairing" and requires the solving of a linear system, as it is usually the case with implicit integration.
Graph Laplacian
A more general way to define the Laplacian operator is to do it on any graph, not only a manifold mesh. And in words, this is a matrix that has the degree (of incidence) of each vertex on the diagonal and then -1 for each connected vertex index and everything else is zero. This can be summarized by the equation
$L = D-A$,
where $D$ is the degree matrix and $A$ the adjacency matrix (1 for each neighboring vertex). But a more interesting way is to use the incidence matrix $B$, which encodes edges in pairs of the form [1 -1]. And this results in a symmetric form factorization:
$L = B B^T$
A mesh is a special case of a graph. So, this tells us that the above Laplacian will always be a symmetric matrix (under normal assumptions) that in turn can come from some kind of "difference" rectangular matrix multiplied with itself transposed. Which makes sens given that tn vectory analysis the Laplacian is the "del" operator or the gradient dotted with itself: $\Delta = \nabla\cdot \nabla$ (sometimes denoted $\nabla^2).
Cloth simulation
Smoothing is equivalent to having a spring on each edge (according to Desbrun et al.), but a special type of spring with zero rest length. This is related to the minimization structure described above. It also explains the shrinking.
We could then build a cloth simulator that reproduces the same behavior as Laplacian smoothing. This could be done using the Baraff-Witkin method (or any other implicit scheme). In this case, the method should be equivalent to "implicit fairing". But it can also be done with an explicit integration step.
Given that a zero-rest length spring doesn't really have a direction, but it's more like a 3 DOF constraint, it is also fully linear. It is equivalent to 3 independent springs acting on each axis. This makes the forces (thus the Laplacian) a linear operator acting on the positions. Thus the stiffness matrix (gradient of force, Hessian of energy) is the Laplacian matrix itself. All the pieces fit!
A first extension that we could think of is to implement a smoothing algorithm based on cloth simulation that doesn't target zero rest length, but only a fraction of the rest length. Of course, it will make things more complicated. For one, the Laplacian operator will no longer be what it used to be, as it's being replaced by the stiffness matrix and nonlinear spring forces. But in my mind it would offer more control over the smoothing process. And the reasoning also goes the other way around: cloth simulation is essentially a mesh smoother. Add on top of that bending springs or volume preservation and you get a powerful geometric tool.
The fact that cloth simulation is a "smoother" is even more clear when doing explicit integration or constraint solving. There is a pattern for accumulating forces in almost any solver, that I call the Jacobi fashion. And this uses the neighborhood of each vertex or incident constraints/springs. You could do it in two passes: one over edges, compute some edge quantities and then sum them up for each vertex (either using the incidence info or not). Or you could do it in one go for each vertex (and its one ring neighborhood), which starts looking a lot like mesh smoothing or a convolution kernel. And in the end a Jacobi constraint solver (or Gauss-Seidel) is just a smoother. In fact relaxation solvers are sometimes called smoothers in the context of multigrid. And there are actual results out there in fluid simulation where Jacobi acts as recursive convolution kernels - they call them Jacobi-Poisson filters.
And to go a bit into details, we already know from the context of implicit cloth simulation that the Laplacian operator corresponds to the stifness matrix $K=k J^T J + K_g$, where $k$ is the springs stiffness (considered the same for all springs), $J$ is the Jacobian of the constraint functions (used by the springs too) and $K_g$ is something that I will call geometric stiffness and ignore here. Here $J$ stands for an equivalent (transposed) of the adjacency matrix $B$ - it's actually an equivalent with extra orientation in space. $J$ also stands for the Jacobian matrix, which is just a generalized and transposed gradient.
Now, if we take a projection operation, which is at the heart of any constraint solver, you get
$p' = (I - J^T (JJ^T)^{-1}J) p$.
This again gives us a glimpse of an equivalent Laplacian operator. But even more, this is equivalent to solving a linear system
$(JJ^T)x = b$
which in fact is what is done in practice for finding out the Lagrange multipliers, here denoted by $x$. And the matrix of this system is clearly a Laplacian of some graph itself and the whole linear system corresponds to a discrete form of a Poisson equation. I posit that this Laplacian is nothing else but the one coming from the dual graph of the constraints/springs. And it is something that I've been glimpsing throughout the years. It's not especially hard if you look closely at a PBD or impulse solver (especially Jacobi ones). But my current realization, which I've been guessing all this time, is that a constraint solver is really a smoothing operation of some field (velocity, position) in the end.
However, in the end, smoothing and cloth simulation are still distinct processes. Maybe they are both more or less equivalent to each other and they both act as spatial filltes. It is also said that the "contangent Laplacian" can be derived from FEM, thus it cannot come from springs alone, but rather triangle elements. So, it's just different flavors of "Laplacians" that appear in mass-spring or constraint-based simulations. Also, smoothing seems to act more on the curvature part of the surface, thus bending forces or constraints may play a bigger role on the side of cloth simulation. In fact, usually they are all super-posed on each other (stretching, shearing, bending) resulting in a global stiffnes matrix or Laplacian. So, there are a lot of unanswered questions in my head and experiments to be done.
Bibliography
Desbrun et al., "Implicit Fairing of Irregular Meshes using Diffusion and Curvature Flow", 1999
No comments:
Post a Comment