vector field

MESH DEVELOPABILITY

Some weeks ago I had an interesting conversation on the Houdini Discord Server about Keenan Cranes fairly new paper Developability of Triangle Meshes. The paper is, of course, excellent as always and the topic quite interesting in general. Nevertheless I somehow missed it when it was officially published and even worse – I completely forgot about it although I was at a conference around two years ago where Keenan Crane gave a beautiful talk about their new method. Fortunately there are people around who are really quick in finding new papers and after it was pointed out to me by Spinynormal I did a quick and dirty implementation in Houdini. It’s just a rather rough prototype without taking too much care about all the details but nonetheless it seems to work quite well.

To be honest, at the time when I saw Keenan Cranes talk I was first a bit sceptical. I really liked the idea and thought wow, that’s brilliant but I also had my doubts about the practicality of their method. Since the result depends largely on the structure of the mesh and not so much on it’s shape (or other properties), flipping some edges, for instance, might lead to a very different solution. And to put it simple, it wasn’t sure if this is an advantage or rather the opposite. Well, now, years later and after playing around a bit with the implementation in Houdini, I’m pretty much convinced that it’s an advantage. It’s quite interesting to watch how the mesh evolves into different shapes by just changing small parts of the triangulation. But then again, I don’t need to use it in practise.

K. Crane, O. Stein, E. Grinspun: Developability of Triangle Meshes

VECTOR FIELDS … AGAIN

I guess I know what you think: “Not again another post about vector fields!” Well, that might be true but vector fields are omnipresent in computer graphics and therefore it’s incredible important to get a basic understanding of vector calculus and how it’s used in Houdini. A vector field is, as the name implies, a field of vectors in n-dimensional space. In Houdini terms, it’s a bunch of vectors assigned to points, vertices, polys or voxels of a mesh or a volume. That’s roughly what a vector field is. It’s easy to understand and easy to visualize. However, much more interesting than what it is, is to know what we can do with it and how we compute it. And that’s exactly what I’m trying to explain in this blog post. I’m not going into all the details but trying to give a general overview in the context of Houdini’s Surface Operators, better known as SOPs.
So, how do we compute a vector field in Houdini? Well, honestly, there is no single answer to this question because it mainly depends on what kind of field we would like to generate. One of the most important vector fields, however, is the gradient vector field. So, what is a gradient vector field? The gradient is the first order derivative of a multivariate function and apart from divergence and curl one of the main differential operators used in vector calculus. In three-dimensional space we typically get it by computing the partial derivatives in x, y and z of a scalar function. In other words, if we apply the gradient operator on a scalar field we get a vector field. Well, that’s more or less the general definition but what does this mean in the context of Houdini? I’ll try to explain it on a simple example.

Lets create a grid, wire it into a mountainSop and we’ll get some kind of a landscape like in the images below. Now let’s have a look at different points on that “landscape” and measure their height. For every of these positions, the gradient is the direction in which the height (P.y) is changing the most. In other words, it is the direction of the steepest ascent. We could check this by simply sampling the height values of our landscape around the points. Just copy small circles onto the points, use a raySop to import the height attribute and search for the point with the largest value. This way we get a good approximation to the gradient, depending on the number of sampling points. Of course, that’s not how the gradient is computed but it might help to get a better understanding of it’s meaning. The height is simply the scalar input for the gradient operator and the result is a vector. If we apply the gradient operator on every point of the geometry, we finally get our gradient vector field over the mesh.

 

So far so good, but what’s so special about the gradient vector field you might ask. Well, a gradient field has several special and important properties and we are going to explore some of them by having a closer look at our previous example. For instance, if we generate contour lines of the scalar field (our height attribute which is the input for the gradient operator) we can see that the vectors are locally perpendicular to the isolines. If we rotate the vectors by 90 degrees, or rather compute the cross product with N, we get a vector field which is tangent to the contour lines. While both vector fields are tangent to the surface,  the rotated vector field is in terms of its properties pretty much the opposite of the gradient. I won’t write much about it since I already did in another post, but nevertheless it’s important to mention what the main difference is. The gradient vector field is curl-free, it’s rotated counterpart, however, is a solenoidal vector field and hence divergence-free. If the field is curl- and divergence-free, it’s a laplacian (harmonic) vector field. But let’s go back to the gradient for now and have again a look at our “landscape” example.

 

Imagine, you’re standing on the top of the hill, put a ball on the ground and give it a light push in a specific direction. Most probably the ball will start rolling and if you’re not fast enough to stop it, it’s rolling down the hillside to the very bottom of the valley. Well, that’s nothing special and not much of a surprise. What is special, however, is the path the ball will take on its way down. It’s the path on which the height decreases locally as fast as possible and this route is exactly following the gradient vector field, or more precisely, the opposite direction in our example. This principle is called gradient descent/ascent and is used extensively for local minimization and maximization problems, as well as for local optimization problems in general. But mind, the keyword here is local. For global optimization problems we have to use different methods. This is quite obvious if we again look at our example. Depending on the sample point, or rather starting point of our path, we reach the local but not necessarily the global minimum/maximum. To make sure, that we actually find the global minimum, in other words, the point at the very bottom of the valley, we have to use a rather large number of samples.
Since the paths are following the gradient vector field, they are perpendicular to the contour lines of our height field and together they form a network of conjugate curves. Conjugate curve networks are quite important for many different applications such as remeshing, polygon planarization, parameterization and so on. The paths following the rotated vector field are, of course, parallel to the contour lines of our scalar function.

 

Talking about the gradient typically means that we have a scalar field, apply the gradient operator and finally get a vector field. The good thing about a gradient vector field is that we can reverse this operation. In other words, if we have a vector field which is a gradient field, we can calculate the original scalar field. This is possible as long as the vector field is curl-free, what a gradient field per definition is.
Understanding the concept of the gradient operator is quite important since it’s interrelated to pretty much everything in (discrete) differential geometry. For instance, if we apply the gradient operator on a scalar function, in other words, we compute the first order partial derivatives we get a vector as we already know. If we compute the second order partial derivative, the result is a matrix called the Hessian matrix. The trace of this matrix is then the Laplacian of the scalar function.

But now back to the beginning of the post. So, how do we compute the gradient? Well, the definition is quite simple: We have to compute the first order partial derivatives in x, y, and z. On volumes this is pretty easy to do since it’s basically a spatial grid and we can translate the equation pretty much one-to-one to VEX and use it in Houdini. On a parameterized surface we have to change it slightly and compute the derivatives in u and v. On meshes, however, Houdini’s prefered geometry type, we have to do it differently. The simplest method is just using the polyFrameSop in attribute gradient style. This way Houdini is doing the work for us and we don’t have to care about the underlying algorithm. If we don’t want to use the polyFrameSop for some reason, it’s fairly easy to implement our own version in VEX. If you take a look at the HIP you can find three examples. The first method uses the fact that the gradient vector field is perpendicular to contour lines. The second method is basically based on a piecewise linear function and the third is using point clouds to compute the gradient.

And finally some examples of using the gradient operator for various applications, such as magnetic field line generation, flow field generation, medial axis extraction and so on.

 

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.

 

 

 

GUIDANCE VECTOR FIELD

A good friend of mine is an artist. Last year when he did some kind of an installation he asked me if I could help him on an interesting problem. For an interactive installation he basically wanted to work with reflections on a three-dimensional surface. The general idea was quite simple and the setup basically as follows:
He wanted to have two sculptures in a room which could change their position in space. Pretty much in the middle of the room there is a three-dimensional surface, which could be deformed by using some servos with Arduino. This surface responds to the movement and position of the visitor and get modified in such a way that the two sculptures are always visible as reflection on the surface. Even if they are not visible directly. Well, that was the plan. He already had great ideas how to build the reflecting surface in reality and also knew how to build the servo-rig. Tracking the visitors wasn’t a problem either but what he didn’t know was how to define the shape of the surface based on the incoming data, which in turn he needed to control the servos. That was the problem I was supposed to take care of since I promised him to help.
So, just to summarize the “surface-problem”:

  • The surface should be deformed in such a way that both sculptures are visible as reflection as much as possible
  • The surface should be as smooth as possible because of aesthetic reasons
  • Due to material constraints, the surface area should stay pretty much the same
  • Since it’s all about a responsive surface in an interactive installation, it has to work fast

So, because I had no clue how to solve this problem I started doing some prototyping in Houdini. I first set up a simple scene with two spheres representing the sculptures, a grid representing the surface and of course, a “visitor” (in fact, another sphere). After some unsuccessful tests I started to compute various vector fields describing visual rays and reflections in the scene. And suddenly it made “click” – I knew how to do it. Looking at the setup I had so far, it became quite apparent that I could just use the vector field as gradient field in the Poisson equation. I tried it immediately and it worked.

Responsive surface
Visual rays

The vector field was easy to compute. For every polygon of the grid I had to compute the vectors pointing to both “sculptures” by using a weight-function based on distance and reflection angle to finally blend between them. This way I got the direction in which visual rays needed to be reflected and hence I also knew the optimal orientation of each polygon. After that I applied a local rotation to each polygon of the mesh and obtained the new orientation. By doing so, every polygon turned to it’s “best” direction and the mesh was torn apart. Subsequently it needed to be reconstructed based on the new orientation of its polygons. And this is exactly the point where the Poisson equation came into play.

The setup was quite simple. Based on the new orientation I calculated a gradient for x, y and z, computed the divergence for each of these and used the result as new gradient field in the Poisson equation. After solving the Poisson equation, which is in fact a sparse linear system, the result is a piecewise continuous function. In other words, it is the final mesh reconstructed in the least squares sense.

 

LINE INTEGRAL CONVOLUTION

Line integral convolution is a texture synthesis method and a very nice way to visualize a vector field. Compared to streamline generation for instance, it has the advantage that no seed points are needed and usually the direction of the field is better readable. Furthermore it can be used to visualize the vector field in the viewport but also for showing it as a render. LIC is a simple but flexible technique and basically works by blurring a noise texture in the direction of a vector field. This is done by using a convolution filter kernel which is oriented along the field, what is fairly easy to do in COPs. The resulting texture is then applied to the UVs of the mesh. It’s

While such texture based methods are most naturally used on two-dimensional vector fields and generally easy to generate, it becomes difficult when the field is defined on a three-dimensional surface. As long as the surface can be parameterized by using two-dimensional coordinates (in other words, mapping the 3D space into a 2D parameter space), it works pretty much the same way. The LIC-texture is generated based on the vector field in parameter space, and finally mapped back onto the original surface. The problem with this method is, that the parameterization typically introduce distortion and the resulting LIC-texture is therefore distorted too. To avoid this problem and to use LIC on any arbitrary three-dimensional surface, the texture is computed in the local space of each triangle. To do this, I’ve simply encoded the vector field and point positions as colors which I used in COPS for the convolution process. The streamline integration, or more precisely, direction evaluation, however, is done in SOPS on the original surface to ensure that no triangle boundaries are visible in the final texture.

The vector field used in the images below is the maximal principal curvature direction. You can find more about it here.

 

TANGENT VECTOR FIELDS ON MESHES

Vector fields on meshes are useful for numerous things and hence they have many applications in computer graphics, such as texture synthesis, anisotropic shading, non-photorealistic rendering, hair grooming, surface parameterization, remeshing and many more. Sometimes they are quite easy to compute and sometimes it gets tricky. Depending on the task, we maybe have to work with an existing vector field but most of the time we need to generate a new one, typically based on a small set of user-defined constraints. Either way, we usually want the vector field to be tangent and more importantly, smoothly varying over the surface. Well, the tangent part is easy. One way how we could generate such a field is simply by computing the gradient of a scalar function, and the result is of course a tangent vector field. This is easy to do but the resulting field might not be what we want. If, for instance, we need to compute a vector field which should rotate around a point, it won’t be possible using the above mentioned method. The reason is quite simple – since it is a gradient field, it is per definition always rotation-free. In other words, curl is always zero. Of course, if we compute the orthogonal complement of each vector, in other words rotate it around 90 degree, we would get rotations. But at the same time we would naturally end up with a divergence-free vector field, which again might not be what we want.

Instead of utilizing the gradient, a much more flexible way to design vector fields is by prescribing user-defined constraints, which are then interpolated over the surface. These constraints are typically sinks, sources and vortices but could also be some edges on the mesh, which guide the direction of the field. Well, at this points things become a bit tricky but fortunately it is a very well researched topic and Google unveils the great wisdom hidden in many research papers. A paper I’ve particularly enjoyed is “Design of Tangent Vector Fields” by Fisher, Schröder, Desbrun and Hoppe. The described method is based on quadratic energy minimization and relies only on standard operators from discrete exterior calculus.

M. Fisher, P. Schröder, M. Desbrun, H. Hoppe: Design of Tangent Vector Fields

After reading the paper I did some prototyping in python and while it worked generally quite well it was of course rather slow. It was probably around 2014 at the time and I wanted to port the algorithm to the HDK or at least use a sparse matrix library in python but something intervened and later I lost somewhat interest in it. Well, now it’s 2017 and this blog seemed to be a good opportunity to dig out the old file. It´s still slow (even on a fast computer) but if you’re not in hurry it does a pretty good job.