# CP MESHES

About a month ago I wrote my first serious post in this blog. It was about circle packings. Well, that’s nothing special because pretty much everybody knows what a circle packing is and many have used it in the field of generative art to produce stunning images and patterns. But what most people don’t know is how many uses such a packing has in the context of discrete differential geometry. Quite interesting for instance is the so called CP mesh. Originally developed to rationalize freeform surfaces in architecture and design, a CP mesh is a triangle mesh whose incircles form a packing.
To compute the CP mesh we first need to set up an optimization problem based on the original mesh. This is done by formulating an energy function measuring the incircle packing property, the distance to the original mesh and the distance to the boundary. After solving the nonlinear minimization problem in the least squares sense we finally obtain new point positions forming the required CP mesh. However, it is essential to bear in mind that it’s not possible to reconstruct a CP mesh from an arbitrary surface. In other words, there is always a deviation to the original mesh no matter what energy function we use. This is clearly visible on the images below. The picture in the middle shows a CP mesh for which the surface proximity is used in the energy function whereas for the picture on the right the boundary proximity is used.

CP meshes have some useful properties which are quite relevant to various algorithms and applications. One of these properties is that they form an ideal primal/dual structure which guarantees orthogonality, non-intersection and convexity. This is important in case we are using discrete differential geometry operators: Orthogonality ensures good numerical qualities of the Laplacian and furthermore it ensures that we can safely compute primal values based on the dual values and dual values based on primal values respectively.

The combinatorial relation between primal and dual mesh is quite important in the context of discrete differential geometry, especially when we are using DEC methods. For tight error bounds it’s therefore essential to work with a good primal triangulation which ensures the above mentioned properties such as non-intersection, orthogonality and convexity of the dual. To achieve orthogonality, we usually place the dual vertices in the circumcenter of its corresponding primal complex. This however could lead to intersections in case our primal mesh contains obtuse triangles. To avoid this problem we could use the barycentric instead of the circumcentre dual. But this in turn would fail to satisfy orthogonality and convexity. We could circumvent this new problem partially by generally using the circumcenter but switching to the barycenter in case the dual vertex is outside of its corresponding primal triangle. What we could do as well, is to project the vertex onto the triangles boundary. However, to cut a long story short, the only way to avoid all these problems is to use a proper triangle primal mesh. And a CP mesh is such a mesh.

(a) Triangle mesh.  (b1, b2) Dual mesh based on the circumcenter of its primal triangulation. Obtuse triangles for which the circumcenter is outside of the triangles boundary are leading to intersections of the dual mesh  (c) Dual mesh based on the barycenter of its primal triangulation. In this case there is no intersection but it fails to satisfy orthogonality.  (d) Dual vertex is projected onto the edge of its corresponding triangle.  (e) Optimized CP mesh  (f) Incircle packing.  (g) Dual mesh based on the incircle packing of the CP mesh

Another important property of a CP mesh is that if we extrude the edges of the dual polygon we get a torsion free structure. In other words, its physical realization could be build out of beams which intersect in a common axis passing through the incircles center. And last but not least, CP meshes can be used to approximate circle packings on a surface. This could be done by either fitting a circle into each vertex gap between the incircle packing or by using the contact points of the sphere packing. Both methods generally produce a fairly good packing. A real circle packing however, is not possible on an arbitrary surface.

A. Schiftner, M. Höbinger, J. Wallner, H. Pottmann: Packing Circles and Spheres on Surfaces

# BALANCING 3D SHAPES

3D printing became rather popular over the last years and many people started to print their own models. The technology is great and the workflow generally quite hassle-free: Just do the design in your favourite 3D software, save the file and then send it to manufacturing.
There is one problem, though (apart from the necessity to have a watertight model): It’s easy to forget about the laws of physics, or more precisely, to forget about gravity. Consequently, the printed model which is intended to stand, simply doesn’t. The reason for this is most probably because it’s unbalanced and thus just tuple down. However, a great solution to this problem is provided in the excellent paper “Make It Stand: Balancing Shapes for 3D Fabrication” published by the igl Zürich. The authors propose a two-step balance optimization algorithm which works by firstly modifying the volume of the object and secondly, by deforming the original shape. Since I like 3D printing and have first-hand experience with the “it doesn’t stand” problem I thought an implementation of the paper might be worth a try.

The method described in the paper basically works as follows:

• Define support polygons of the shape
• Calculate center of mass
• Carve out inner void based on computed center point
• Calculate new center point and check if its projection along gravity is on support polygons
• If center point isn’t on support polygons, furthermore deform original shape to ensure stability

Computing the center of mass (not to be confused with the centroid) is easy. It can be done quite efficiently by using the divergence theorem. It’s described in the paper but in case you are interested in details, check out the link. For the mass carving part Houdini already has the perfect toolset, namely VDBs. After converting the mesh to a VDB I used vdbReshapeSdf to compute the overall offset of the surface. The final inner volume is then computed iteratively during optimization in a volumeWrangle based on a mask. For the shape deformation I used on an old custom SOP based on laplace editing. And finally, the overall minimization problem was solved in Scipy within a SopSolver. Even though this setup was quite inefficient since all matrices have to be rebuild in every iteration, it worked surprisingly well. Of course, there is much room for improvements but tbh. I was more interested in the general method than in the exact implementation of the paper.

R. Prèvost, E. Whiting, S. Lefebvre, O. Sorkine-Hornung: Make It Stand: Balancing Shapes for 3D Fabrication

# 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

# MINIMAL SURFACES

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. For the examples below I had no problems to interactively modify an average dense mesh in realtime on my rather 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.

U. Pinkall, K. Polthier: Computing Discrete Minimal Surfaces and Their Conjugates

# SURFACE CURVATURE

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:

G. Taubin: Estimating the Tensor of Curvature of a Surface from a Polyhedral Approximation
D. Panozzo, E. Puppo, L. Rocca: Efficient Multi-scale Curvature and Crease Estimation
M. Mayer, M. Desbrun, P. Schöder, A. H. Barr: Discrete Differential-Geometry Operators for Triangulated 2-Manifolds

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.