If you ever had to compute an object oriented bounding box by yourself you probably came across the term covariance matrix. Well, this is no accident since covariance matrices are widely used in computer graphics for many different applications such as data fitting, rigid transformation, point cloud smoothing and many more. A covariance matrix is a N x N matrix basically measuring the variance of data in N dimensions. In case of Houdini’s three-dimensional coordinate system it’s therefore a symmetric 3 x 3 matrix capturing the variance in its diagonal and the covariance in its off-diagonal. Computing the covariance matrix is fairly easy and could be done efficiently in VEX or Python/Numpy. The interesting thing about the matrix is that it represents the variance and the covariance of the data. In other words, it describes the “shape” of the data. This means that we can, for example, find the direction of the largest variance measured from the centroid, which consequently is the largest dimension of the data. This is done by decomposing the matrix into its eigenvalues and eigenvectors. In case we are working with a 3 x 3 covariance matrix we get three eigenvectors. The eigenvector with the largest eigenvalue is the vector pointing into the direction of the largest variance. The eigenvector with the smallest eigenvalue is orthogonal to the largest eigenvector and, of course, pointing into the direction of the smallest variance.


Finding the eigenvectors for a small matrix is quite easy and could be done, for example, by using singular value decomposition (SVD) in Numpy. In VEX we don’t have such a function but we could implement a (inverse) power iteration algorithm to find the largest and/or smallest eigenvector. The missing third eigenvector is then computed simply by talking the cross product of the other two already known vectors.
Finding these eigenvectors is quite useful for many applications. For instance, if we need to compute object oriented bounding boxes for each polygon of a mesh. We could do this by building the covariance matrix for the points of each polygon and then searching for the largest eigenvector. Together with the normal, which is the smallest eigenvector, we finally can compute the oriented bounding box. Another useful application of the covariance matrix is the computation of point normals for a point cloud. In this case we need to find the smallest eigenvector, which represents the best fitting plane and thus the normal of the point. Building the covariance matrix is also important if we want to find the rigid transformation from one object to another. All these are just a few examples and to cut a long story short, covariance matrices are quite useful and definitely worth looking into.





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


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


Non-photorealistic rendering is something I have been quite fascinated with for a long time. Unfortunately I never had enough time to really get into it. But when a good old friend asked me if I could help him on one of his projects I thought this would be a good opportunity to finally dig deeper into NPR land.
My friend is one of the few persons who pretty much never uses a computer. And for his work he doesn’t have to. He is a brilliant artist making wonderful drawings mainly with colored pencils. Around six to seven years ago, he started one of his newer projects. He developed his own alphabet as a kind of abstract notation which is based on rectangles and four different colors. In the following years he made hundreds (or maybe even thousands) of drawings which show all kinds of variations of this notation. Usually these drawings are translation of one or two words into his alphabet. At some point he finally wanted to translate more than a few words but since he is doing all by hand, this became rather time consuming. One evening, when we sat together over a drink we talked about his project and he told me about his idea of translating “The Pacific”, a chapter from the famous book Moby Dick into one of his drawings. It wasn’t that much text but since he also wanted to try out different variations it became rather tedious to draw it by hand. In Houdini however, this could be done in a fraction of a second. To test it out I quickly threw together some nodes and with the help of python and the mighty copySop I had a working setup shortly after. After that I simply used Wren for rendering out the wireframe as Postscript file which could then be plotted on paper.

My friend was quite happy with the result. It did save him a lot of work and time and now he could concentrate on coloring and finalizing the drawing. I for one, however, was asking myself if it was even necessary to hand colorize the prints. Wouldn’t it be possible to just write a shader and let Houdini and Mantra do all the work? I was rather confident and so I started working on the shader. On one hand it was important to simulate my friends drawing style as close as possible and on the other hand the shader had to be as flexible as possible. I first thought I could scan parts of his drawings and simply use it as texture but this wasn’t neither flexible enough nor did it look natural. It was quite the opposite, it looked terribly unnatural. After some more tests I finally decided to use only layered procedural textures. This made rendering rather slow but since I had to render just a few large stills, this wasn’t much of a problem. At the end it worked quite nicely and using fully procedural textures made it possible to simulate many different drawing styles.