Computing curve skeletons of 3D shapes is an important tool in computer graphics with many different applications for shape analysis, shape matching, character skeleton construction, rigging and many more. It’s an active field of research and over the years, numerous algorithms have been developed. Most of these algorithms are based on mesh contraction using smoothing or thinning. Another popular approach is to use the voronoi diagram generated by the points of the mesh. Quite an interesting alternative, however, is to rely on spectral analysis for computing the skeleton. The idea is fairly simple and works as follows:

Compute the first non-zero eigenvector of the Laplace-Beltrami operator, namely the Fiedler vector.

Some time ago I wrote a blog post about the covariance matrix and how it’s eigenvectors can be used to compute principal axes of any geometry. This time I’m writing again about eigendecompositions but instead of the covariance matrix it’s the Laplace-Beltrami operator which is decomposed into its eigenvalues and eigenvectors. Now, the question is why would we like to find eigenvalues and eigenvectors of the laplacian? Well, there are a number of applications for which this is very useful. For instance, calculating the biharmonic distance on a mesh. Dimensionality reduction and compression are obvious examples but beside that it’s quite important for various other things too, such as parameterization, smoothing, correspondence analysis, shape analysis, shape matching and many more.

By projecting the point coordinates onto the eigenvectors we can transform the geometry to the corresponding eigenspace. Similar to the Discrete Fourier Transform, the new mapping could then be used for doing various calculations. which are easier to perform in the spectral domain, such as clustering and segmentation.

Eigenvector 3

Eigenvector 6

Eigenvector 8

Eigenvector 10

Eigenvector 15

Eigenvector 30

Eigenvector 50

Eigenvector 1-100

A very important property of eigenvalues and eigenvectors of the Laplace-Beltrami operator is that they are invariant under isometric transformations/deformations. In other words, for two isometric surfaces the eigenvalues and eigenvectors coincide and consequently their eigenvector projection also coincide. Since a character in different poses is usually a near-isometric deformation, eigenfunctions are quite useful for finding non-rigid transformations or for matching different poses. For instance, if we have two characters like in the example below, we can rely on eigenvectors to compute a correspondence matrix and therefore find corresponding points on the two meshes. Subsequently, these points could then be used to match one pose to the other. Now you might ask why do we need eigenvectors for this? Well, if it is the same character in different poses, or more precisely, the same geometry with the same pointcount, this is easy. If it’s not the same geometry, things get much harder. This is typically the case if we are working with 3D scanned models or shapes which have been remeshed like in the example below. The two characters have different pointcounts: The mesh of the left character has 1749 points, the mesh of the right character has 1928 points.

A very useful tool is the Fiedler vector which is the smallest non-zero eigenvector. It is a handy shape descriptor and can be used to compute shape aware tangent vector fields, to align shapes, to compute uvs, for remeshing and many other things. The images below show that isolines based on the Fiedler vector are nicely following the shape of different meshes.

When I saw Entagmas great tutorial about minimal surfaces I remembered an old file I did to help a friend some years ago. He wanted to generate a smooth surface between some given boundaries and it should behave similar to elastic textile. To simulate it he used the same method which had been invented about 40 years ago by Frei Otto. He built a physical model out of thin wire and dripped it into soap film. After taking it out of the fluid the result was a perfectly smooth surface in the form of a soap film spanning between the boundaries.
After some experimentation and many models later, the question was finally how to digitize the surface. 3D-scanning was obviously not possible. He tried to model it based on photos but that didn’t work out very well. Repeatedly subdividing a coarse mesh wasn’t successful either. Thereby the solution is quite simple. Don’t try to digitize the surface, just digitize the whole process and skip the physical model.
Mathematically such a surface is a minimal surface or more precisely, it is a surface of minimal area for given boundary conditions. Another important definition is that it has zero mean curvature at all points. What he had to do is to minimize the overall area of the surface until the mesh reaches equilibrium. The resulting surface is a minimal surface, or at least a very good approximation of it. So, what’s the best method to compute such a surface?
Well, one way is of course to use a spring model, set the rest length to 0 and run the simulation. That’s basically what Entagma did in their tutorial. The problem with this method is that it´s rather slow and works best only on a fairly regular and clean mesh. On irregular meshes it’s highly likely that it won’t converge to a minimal surface.
A probably better method is laplacian smoothing. Again this works best on a regular mesh but it’s readily adaptable by using proper weights. While both methods might be slow due to the fact that they use an iterative approach, it’s generally better to rely on laplacian smoothing. This will result in a better approximation of a minimal surface.
Of course, there is a third method which again relies on the laplacian but uses a slightly different approach to solve the same problem. It´s based on the paper “Computing Discrete Minimal Surfaces and Their Conjugates” by Pinkall and Polthier and works by minimizing the Dirichlet energy. This leads to an efficient numerical method to compute the desired minimal surface. The basic idea is simple – just move the vertices of the mesh in a specific direction for a specific distance and the result is a minimal surface.
To do this, we first need to find the direction in which the area changes the most. In other words, we need the gradient of the area. Apparently this is the mean curvature normal which we get fairly easy by using the Laplace-Beltrami operator. The length of the mean curvature normal is the mean curvature and as we already know, it should be 0 for minimal surfaces. With this information it´s easy to set up a linear system in which we set the curvature to 0 and simply solve the system of equations for the surface with harmonic coordinate functions. The solution we get are the new point positions for the minimal surface we were looking for, or at least a very good approximation to it.
When using a fast sparse matrix solver this method is not just quite accurate but also fast. I’ve used Eigen as linear algebra library in the HDK for this purpose and had no problems to interactively modify an average dense mesh in realtime on my 5 year old laptop.

The same setup could be used for membrane and tension surfaces since they share many similarities with minimal surfaces, which is quite nice for playing around with the shape.

Curvature is a fundamental tool in computer graphics and important for many techniques such as remeshing, smoothing, simplification, npr-rendering and many more. Unfortunately Houdini doesn’t provide something like a CurvatureSop out of the box. Of course, there is the MeasureSop which can calculate some kind of curvature but the result is far from being optimal. It doesn’t differentiate between different types of curvature, such as gaussian or mean curvature nor does it give you the important information if it’s positive or negative. In other words, if the shape is concave or convex.

Since I needed to calculate curvature quite often in the past I decided to build my own CurvatureSop. Well, one might think that this should be quite easy but in fact it isn´t. The problem is that in Houdini we are mostly working with meshes. Meshes are nice to work with for all sorts of things but they are not that great if we need to compute derivatives. And computing derivatives, or more precisely, the second derivative is exactly what curvature is about. In the continuous case this is simple and straightforward to do as it is pretty much on volumes (structured spatial grids). On an arbitrary mesh however, we face the problem that the second derivative of a piecewise continuous function is pretty much everywhere null, at least if we are working with triangulations. The same is not necessarily true on quadrilateral meshes but to keep it simple and not becoming even worse, we´ll stick with triangles. On a triangular mesh, curvature exist at vertices or along edges and therefore it’s highly dependent on the quality of the mesh. One way to partially minimize this problem is to locally fit a polynomial to the points around the vertex and use the fitted surface to compute the curvature analytically. While curvature estimation using this method is quite accurate, it is also very sensitive to the distribution of points in the neighbourhood and thus still dependent on a regular triangulation. A probably more straight forward way to measure curvature is to apply the Laplace Beltrami operator directly on the triangular mesh. This works generally quite robust as long as we are working with non-obtuse triangles. To cut a long story short, there exist many different algorithms to compute curvature on meshes but all have their pros and cons very much depending on the mesh. To compensate for this I finally implemented three different methods and usually use what’s giving me the best result for the mesh I´m working with. All three algorithms are using a different approach and are based on the following papers:

Beside using different algorithms, the CurvatureSop allows to compute the main curvature types, such as mean curvature, gaussian curvature and principal curvatures as well as principal curvature directions.