# THE JOY OF USING MATRICES

## Fom VEX/OpenCl to sparse matrices

In my last blog post I wrote a bit about the possibility to (mostly) circumvent the need of using sparse matrices in Houdini and instead setting up and solving the system in VEX/OpenCl.  This time, however, it’s about the other way around and how a typically VEX/OpenCl based setup could be implemented and solved in matrix form.

Pretty much everybody using Houdini has at some point done basic Laplacian smoothing on a mesh (with uniform weights), perhaps without even knowing it. It’s a widely used algorithm and fairly easy  to implement by averaging the values of the 1-ring around the current point. In fact, in Houdini it’s even easier than that and we can just use the attribBlurSop. AttribBlur is using OpenCl under the hood and relies on the basic graph Laplacian to blur/smooth point data, such as P for instance. This is done by performing the smoothing/averaging step a number of times which is conveniently exposed as iteration slider in the UI of the attribBlurSop.

However,  this is of course not the only way how Laplacian smoothing could be implemented. Instead of iteratively applying the smoothing step per point we could also set up a simple linear system and solve it subsequently for the final point positions. This is what I’ve done in the file below. It’s a simple example where I’m using a direct solver in Scipy which is basically solving all the equations simultaneously in one go.

In practice it doesn’t make much of a difference which method we use because both will lead to the same result in pretty much the same amount of time (at least for our example). But since Laplacian smoothing is already provided to us by the attribBlurSop it’s most probably best to just use what’s already there instead of implementing it in matrix form.

Anyways, apart from practicality I hope this simple example, together with my last post, helps to clear up some misunderstandings and the question if it’s absolutely necessary to use sparse matrix libraries in Houdini or not.

FILE

ps.: To run the example file without problems you’ll need to have Scipy installed on your system/computer.

# FIEDLER VECTOR AND THE CURVE SKELETON

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.
• Compute contour lines based on the Fiedler vector
• Find the center point for every contour
• Fit curves to the points

# SPECTRAL ANALYSIS

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.

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.

H. Zhang, O. Kaick, R. Dyer: Spectral Mesh Processing
V. Jain, H. Zhang: Robust 3D Shape Correspondence in the Spectral Domain
S. Dong, P. T. Bremer, M. Garland, V. Pascucci, J. C. Hart: Spectral Surface Quadrangulation
B. Vallet, B. Lèvy: Spectral Geometry Processing with Manifold Harmonics
A. Shtern, R. Kimmel: Matching the LBO Eigenspace of Non-Rigid Shapes via High Order Statistics

# VECTOR FIELD DECOMPOSITION

Some time ago I wrote a blog post about generating tangent and smoothly varying vector fields  on arbitrary meshes.  The method I’ve used works basically by specifying some constraints such as vortices, sources, sinks and guiding curves/edges, which are then interpolated over the surface to compute the final field. Depending on the constraints we achieve very different kinds of vector fields which have different properties and hence different applications. So, while the resulting vector field might be exactly what you need as velocity field in a fluid simulation, it might be completely useless for remeshing. The main reason for this is related to curl and/or divergence. Generally speaking, whenever we are computing or just working with vector fields, we always have to deal with at least one of the three interrelated main differential operators – gradient, divergence and curl.

For instance, applying the gradient operator on a scalar field results in a curl-free vector field. If we then apply a local rotation around 90 degrees, we get the orthogonal counterpart in the form of a divergence-free field. This is e.g. exactly what we need as velocity field for a fluid simulation since we usually want fluids to be incompressible. Furthermore we have to take care of open edges at which we usually want to have tangential boundary conditions, at least in case we are working with fluids.

Well, though all this sounds easy, the problem is that in most cases we just don’t have these nice velocity fields. Normally they are neither gradient fields nor solenoidal fields but rather a mixture of both. In other words we are dealing with vector fields which aren’t completely divergence-free. That’s the reason why we need to apply pressure projection in DOPS which finally computes a proper velocity field needed for fluids. This works nicely in the simulation environment on volumes. But what if we need to do the same in SOPS and on surfaces? Of course we could still rely on volumes, do all the operations we need to do and bring the result back onto the original mesh. While this works to some degree, it is far from being optimal – the result will most probably be inaccurate and we will be limited to solenoidal and gradient vector fields. What we could do instead is to work directly with vector fields on the surface. Typically any vector field on a simply-connected domain could be decomposed into the sum of an irrotational (curl-free), a solenoidal (divergence-free) and a harmonic (divergence-free and curl-free) field. This technique is known as Hodge-Helmholtz decomposition and is basically achieved by minimizing the energy functionals for the irrotational and the solenoidal component of the field by solving the Poisson equation. Finally the harmonic component will be computed as the residual.

While divergence-free vector fields are essential for any kind of fluid simulation, curl-free and harmonic fields are essential for many other applications, such as remeshing, texture synthesis and parameterization to name a few. The irrotational component of the vector field is curl-free and hence the field is integrable. It is equivalent to the gradient of a scalar potential function on simply-connected domains. This means that we can reproduce the scalar field from the irrotational vector field because it is exactly the gradient of the desired scalar field. The isolines of the scalar field are therefore exactly orthogonal to the vector field. In case the vector field is not curl-free, we can treat it as a gradient of a function and minimize the difference between the gradient and the given vector field in the least squares sense. Again this leads to the well known Poisson equation.

# AS-RIGID-AS-POSSIBLE SURFACE EDITING

In Houdini we are mainly working with meshes and quite often these meshes need to be deformed. For this reason we can use several tools such as cages or proxies together with the LatticeSop or the newer PointDeformSop. However, sometimes it’s handy to deform the shape directly. Of course, there is the EditSop but getting smooth deformations can be tricky. Even if soft radius is activated, you might end up with unsatisfying results, especially when rotations are involved. This might be the case because the EditSop uses just a spherical distance falloff and no geometry-aware distance, or maybe because local details don’t get maintained during editing, or maybe it’s because the deformed shape looses to much local volume/area compared to the original. Anyway, one of the best techniques to avoid pretty much all of these problems is “as-rigid-as-possible” deformations and it is described in detail in the paper: “As-Rigid-As-Possible Surface Modeling” by Olga Sorkine and Marc Alexa.

A perfect mesh editing tool should enable us to select some vertices to move and rotate them around. The deformation of the handle should get propagated smoothly to the rest of the geometry. And of course, this should happen without removing or distorting surface details and obviously as fast as possible. In other words, the deformation should be smooth and properties such as area, curvature, distance and so on, should be preserved as much as possible. Well, the problem is that if local mesh properties should be preserved, only rotation and translation should be involved since scaling would destroy the local structure. This however is clearly not possible in the case of mesh editing because as soon as we select a vertex and transform it, some parts of the mesh will get scaled and/or shifted. As consequence we´ll end up with distorted surface details.  “As-rigid-as-possible” however, presents a nice solution to this problem. It tries to preserve local rigidity as much as possible. This means, the shape is locally preserved and smooth as well.

Basically this is done in a two-stage approach. First, getting an initial deformation guess by propagating the handle transformation to the mesh. And second, finding the best deformation to maximize the local rigidity by minimizing an energy function. The energy function is based on the one-ring neighbourhood around vertices and measures the deviation to the rigid transformation of the cell. By summing up the energy functions for all vertices/cells of the mesh we get the global energy which is to be minimized to get the final deformation.

The number of iterations required to retain local properties depends mainly on the size of the mesh. Generally, just a few iterations (10-20) are needed until convergence.

The implementation is done with the HDK but instead of using Houdini’s own matrix solver I´ve used the linear algebra library Eigen. For more information about “as-rigid-as-possible” check out the following paper:

O. Sorkine, M. Alexa: As-Rigid-As-Possible Surface Modeling

# BIHARMONIC DISTANCE COMPUTATION

Some time ago I wrote a blog post about geodesics and different methods how they could be computed in Houdini. Geodesic distances are great since they represent the distance, or more precisely, the shortest distance between points on a surface. While this is an important property and useful for many things we might run into a situation in which other properties are even more important, for instance smoothness. This might be the case if we want to generate a smooth vector field based on the distance or if we are working on meshes with holes for example. However, a nice solution to these problems is to compute biharmonic distances instead of geodesic distances. Biharmonic distance is not necessarily the true distance between points on a surface but it is smooth and globally shape-aware.

Contrary to geodesic distance, the biharmonic distance is not computed directly on the surface but on the eigenvectors of the laplacian. Basically it’s the average distance between pairwise points projected onto incrementally increasing eigenvectors. The general accuracy is therefore highly dependent on the number of computed eigenvectors – the more the better. If we need to get exact distances we have to compute all eigenvectors or in other words, we need to perform the computationally rather expensive full eigendecomposition of the laplacian. For a sufficient approximation however, we usually can get away with the first 20-100 eigenvectors depending on mesh density and shape. For the implementation in Houdini I’m using Spectra and Eigen together with the HDK. Even though it’s quite fast to perform the eigenvalue/eigenvector decomposition, the computation of geodesic distances is of course much faster, especially on large meshes.

Y. Lipman, R. M. Rustamov, T. A. Funkhouser: Biharmonic Distance

# GEODESICS ON MESHES

Measuring geodesic distances and the computation of paths on a mesh is widely used in computational geometry for many different applications ranging from mesh parameterization and remeshing to skinning and mesh deformation. After a longer conversation and chat about this topic, and since i thought it could be useful at some point I finally got my feet wet on different algorithms to compute such geodesics in Houdini.
The first and obvious choice was to compute the shortest path between the source point and all the other points on the geometry. This was easy and fast to do but unfortunately it doesn’t compute geodesics. Since the underlying algorithm of the ShortestPathSop (Dijkstra) relies on the edges of the mesh, it computes just a rough approximation of the real distances. Of course, the inaccuracy could be minimized by iteratively smoothing the paths but still, the result will be far from being optimal and the method slow as well.

To circumvent this problem, I tried another method. Instead of measuring the length of the shortest paths directly it relied on the fact that a real geodesic path unfolds into a straight line. So, what I basically did was:

• Compute the shortest path between point A and point B
• Unfold the path onto a plane
• Measure the Euclidian distance between point A and point B on the unfolded path.

This method worked quite nicely but it has one big drawback – it doesn’t work on meshes with holes. As soon as the line drawn between the start- and endpoint of the unfolded path crosses a hole, the result is no longer the true geodesic distance.

The third algorithm I´ve used was the well known Fast Marching Front Propagation Method. You start at the single source point and compute the distance to the 1-ring neighbouring points around it. By doing this iteratively and summing up the distances you´ll finally get the geodesic distance to all points in the mesh. This works because it is easy to calculate the distance to the third point of a triangle as soon as the other two distances are known.
For the first quick implementation i used python but it was way to slow to be usable. After replacing the slow list lookups in the code i got a 3 times speedup but it was still to slow which was the reason for porting the algorithm to the HDK. After that, it worked quite fast (around 270 times faster) but the general problem with this algorithm is, that it doesn´t work with obtuse triangles. In other words, it needs to work on a fairly clean and regular triangulated mesh. If this is not the case, you´ll end up with errors in the distance calculation which leads to larger errors and so forth. To prevent this from happening I had to insert conditional statements for obtuse triangles as well for all other special cases and hence the code got much longer and somewhat cluttered. Usually, when fast marching front propagation is used on 2D and 3D grids, this is not a problem but for arbitrarily triangulated meshes it is.

Finally I decided to implement the great paper: “Geodesics in Heat” by Keenan Crane. The method described in the paper works basically as follows:

• Integrate the heat flow for a source point across the mesh
• Compute the normalized gradient for the heat
• Solve the Poisson equation to get the actual distance

One of the great advantages of the heat method is its brilliant simplicity. Basically it relies only on solving a sparse linear system what could be done quite efficiently by various linear algebra libraries (I’ve used Eigen for this). This makes it possible to compute the geodesic distances from more than one point in pretty much the same computation time, quite contrary to the marching method for example.