Voronoi diagrams have a long history with computer graphics. For me personally, the encounter that first taught me what Voronoi could be used for and more importantly, how to mispronounce the word, was through the Cinema4d plugin, Thrausi.
Around 8 years ago, having just freshly torrented the Cinema 4d* and wanting to copy some sweet internet tutorials I had found, I set off towards fracturing some text with Thrausi’s Voronoi Fracture feature. First though, I had to have a real heart to heart with my parents; I had never illegally downloaded a program before Cinema 4d and freaked the fuck out. Legit I thought that SWAT + Adobe + Maxon would burst through my parent’s front door. My mental ultimatum was either tell them, or something greater would punish me.
Wrong choice. With tears shed, a bruised posterior and a freshly cleared conscience, I booted up some sweet YouTube tutorials. No where else can one find more broken English paired with such high quality knowledge bombs. However Voronoi quickly showed up in all related videos; I was obsessed.
So why was I obsessed? Because I’m a freak. But also Voronoi diagrams look incredibly neat, and not only that, they’re also very good at partitioning data, and by extension, point clouds and meshes.
Two important things to note before moving onto the actual goodness that is Voronoi diagrams. First my previous blog post on Dual Graphs is required reading, if you want to fully understand this post. The second thing I want to make it clear, is that this post is definitely going to go super in depth on Voronoi stuff. And because of that, I highly recommend checking out the Houdini Wiki to brush up on your vex skills, before we get to far: SOAK IT ALL UP
I hope this isn’t too technical too quick, but the way I think about what a Voronoi diagram even is, is first take a surface S (in our case a grid in Houdini), and some points scattered on that surface, P, with a random color on each point. Find the nearest colored point in the scatter P, to each of the points making up our grid S. Finally color S based off the nearest P color.
To abstract further, that basically means a Voronoi diagram is a method of partitioning anything, in this case point positions, by the distance to another point in a separate data set. But it doesn’t need to be positions. As long as you can come up with a distance function for your given problem, you can make a Voronoi diagram. So really one could make a Voronoi diagram all in excel with no visuals. It’d suck and wouldn’t get dem precious Medium reads (FEED ME ATTENTION), but you could definitely do it.
We use Houdini though, so visuals are par for the course. Plus we have some incredible tools at our disposal for finding near points, specifically the near points tools and the point cloud tools.
Enough waiting! Let’s make some goddamn magic.
Here’s a link to a project file (CLICK ME) that we’ll be building off of in each example, for this pass let’s dive into the geometry object “Example 1.”
I’m ready for the hardest coding challenge of my life.
int np = nearpoint(1, @P);
@Cd = rand(np);Wait what.
That’s it?
That’s how easy it is to get an ultra basic Voronoi cell pattern in Houdini. Two lines of code.
Wow! We can check Voronoi construction off our list, maybe add it to our resume. But Jake, this looks like Boring-Town, USA. I made a sad face with it, because it’s so boring :(
Not only that but the output from this example is color on a surface and not actual geometry that can be used for explosions or other goodness. Luckily that’s where that obsession I mentioned earlier comes in.
Remember how I said Voronoi diagrams are the big brother of a dual graph, well from this example, that probably makes no sense. However, I’ve developed a really neat method that *easily* (that’s a lie) shows the deep connection between the two concepts, while building off of this ultra simple method of constructing diagrams using magical colors.
So the first thing to understand about the connection between the two concepts, is that the dual, with a little bit of spiciness, of a very VERY specific type of triangle mesh, is your Voronoi diagram. The tessellation in question is what’s known as a Delaunay triangulation. Delaunay triangulations can basically be defined as a triangle mesh, where when you create a circumcircle around any given triangle, no other points in the mesh, will fall inside that circle. Since the concept of a circumcircle does not come up literally ever in daily life, I made this neat lil animation to help explain what it is, and how it’s affected by the size and shape of a given triangle.
Kinda boring, but it gets the point across. All points of the triangle, must live on the edge of the given circle for the circle to be considered circumscribing the triangle.
Notice from the motion of the centers that the circumcircle center’s movement takes a much smoother path based off the motion of the animation, compared to the triangle center.
At this point it might not be instantly apparent as to why we need a Delaunay Triangulation. Well take a look at what happens when I force the turn the triangle into a sliver. Notice how the bounds of the circle tend off towards infinity? As your triangle gets thinner and the 3 points get closer and closer to living as a line, your circle needs to stretch infinitely far to contain them because thinking back to the middle school geometry days tells us that:
A). That we were a bad student in middle school.
B). That a circle whose radius is infinite, is a flat line.
So the only way to circumscribe a line, is by making a circle large enough (by large enough I mean infinitely big) to create a flat line to contain the points with. To be completely honest, it’s one of those things that definitely needs to be be seen first hand, so in the example file, check out the little “Circumscriber” geo container to play with this idea.
And because the circle is stretching to infinity, that means the radius is infinite; therefore the center is infinitely far from our triangle, which is not going to connect well with circumcenters that don’t tend off towards infinity either. That’s why in the case of a Delaunay triangulation, it’s purposefully built to avoid slivers by ensuring that none of the other triangle’s points can live inside any given circumscribed circle.
So then this is easy:
Step 1). Acquire a mesh that’s been Delaunay Triangulated.
Step 2). Compute the dual.
Well not completely, that’d give you a regular dual graph (that looks pretty similar to a voronoi diagram), assuming we use the method of half edges that we defined last blog post.
The extra spiciness I mentioned above, is the fact that rather than connecting to our neighboring triangle’s centroid, we’re going to want to connect to the center, of the triangle’s circumcircle. That right there is the most important part, and is what gives Voronoi graphs such cool motion. If we scroll back up and check out the motion path animation, we’ll see that the center of the circle (defined in yellow), takes a far smoother path than the center of the triangle (defined in red). So connecting through those will give us a far more robust and stable platform for animated cellular meshes.
The last thing to mention, is that the dual of a Voronoi diagram, is a Delaunay triangulation. That’s what leads to a lot of confusion, because the dual of a Delaunay triangulation isn’t always a Voronoi diagram, remember that requires computing the dual on the circumcenters.
Now if we really wanted, I could stop here and run through the circumcenter formula using the file from last post, along with the beautiful Triangulate2D SOP from Houdini. That SOP when given a scatter of points creates a Delaunay triangulation of that given scatter. However as the name suggests, it’s a 2D operation, and I’m a 3D wizard, so that’s lame as hell. We aren’t going to do that, and also there’s no challenge in that since we’re already like 99% of the way there.
If you want to take the easy way out check out this ODForce thread: http://forums.odforce.net/topic/29958-computing-the-dual-with-vex/#comment-168535
Let’s instead build our own Triangulate 3D SOP, using the pretty color near point method we described above. The last thing I had mentioned in my description of a Voronoi diagram, is that the dual of it, is a Delaunay triangulation. So it stands to reason that we could “compute the dual” on our magical color method, to get an approximate Delaunay Triangulation, from our approximate Voronoi diagram.
The idea behind the system we are about to build is pretty simple, just like in dual’s blog, we’re going to need a way of testing neighbors, however instead of testing neighboring polygons, we want neighboring Voronoi cells. Once we’ve found the neighboring cells, we basically compute the dual on them, to give us a Delaunay tessellation.
Another way to think about neighboring cells, is that your neighbor cell is going to be the second closest test point, so if our voronoi cell in a nearpoints array is index 0, then our neighboring cell point, should in theory, be index 1.
Now that in theory works, and assuming your input object has enough points to sample with, there wouldn’t be any problems with this method. So for shaders it’s a perfectly valid method (see Inigo Quilez’s blog on voronoi edges), but since we’re thinking about meshes, it’s not a foolproof solution. We’re going to need to be a bit more clever about how we find our neighboring cells due to error that can arise from an a solution like this.
Triangulate3D
Let’s hop into the second example container, and get it into this shit. I wrote some psuedocode we’ll be using for this method to help aid the process of implementation. If we’re already at wizardry levels as a technical artist, i’d recommend stopping here and implementing the pseudo code on your own, without following my instruction, as it’s always good practice:
I know it’s not typical psuedocode, but I prefer to write psuedocode as if I were talking it out. Especially since Houdini runs in parallel (see matt’s wiki), so it can be hard to write psuedocode for, since you don’t need to loop through every point manually for each of these operations. A point wrangle will do that for you.
The entirety of this example, until the meshing process, will live inside of one of the For Each Blocks in feedback mode, with “By Count” set as the iteration method. For testing purposes make sure the iterations are set to one, until you’ve completed the walk through and are ready to test. Otherwise we’re going to make this pretty hard on ourselves…
Our first wrangle in the example is going to contain the simple Voronoi generation method we covered initially. However rather than writing to our color variable, we’re just going to want to store the closest point’s ID.
int np = nearpoint(1, @P);
@cell = point(1, "id", np);The last part of the first step would be a way to get our list back from a previous iteration.
Looking at our psuedocode we’ll notice the whole operation returns the vPoints, or the points we’re sampling to construct the diagram with, as opposed to returning the base pointcloud/mesh. And judging from the psuedocode as well, we’re going to be iteratively adding to a list of neighbors. So we need to grab the list off of the same point sampled by nearpoint(), or the points that are feeding into our second wrangle input.
int np = nearpoint(1, @P);
@cell = point(1, "id", np);
int nb[] = point(1,"nearp", np);if(len(nb) > 0){
i[]@nearp = nb;
}
Sick, our first wrangle is all set, and we’re now ready to move onto this portion of our psuedocode:
From each point, A-p1, in any given cell, find the nearest point A-p2 that’s contained in any other cell.
Store the current cell of A-p2 on point A-p1. This is our neighbor cell
This is where we diverge from just using second nearest vPoints point. Rather than just finding the second nearest point, we’re actually just going to test all the surrounding points, in our initial scatter after creating our Voronoi diagram. That way we minimize any potential error that can arise from not having enough surface samples.
So what does something like that look like? Assuming we’re familiar with Houdini, our gut instinct might be to use a For Each Block. However, that can be incredibly slow since in doing that we’ll need to create two geometry streams, one for the current cell, and one for the rest of the points we’re comparing against. Luckily, Houdini is on some next level shit, and they clearly thought about this problem before.
With a function like nearpoints() or pcfind(), if we read the docs, we’ll notice some variations on the function call for a ptgroup (point group). What’s awesome about that, is that feature allows for point group expressions as well. If we think about our cell attribute as point group, it should be pretty straight forward to select the inverse of that. Unfortunately nothing is ever that simple, but with the the help of the great David Chow, we were able to come up with the simple string needed to search through all the points not in our cell.
The key to using the ptgroup expressions, is to store the group as a string variable before using it in the function, otherwise you’ll experience a whole host of strange bugs. Here’s what the string is going to look like
int temp = @cell;
string group = “\@cell!=” + itoa(temp);So to break that down, we’re writing our current cell to a temporary variable, and then in our string we’re converting our cell attribute into a string that Houdini can read, by using a backslash (also known as an escape character) to allow Houdini to read the @ symbol in a string combined with the standard != does not equal operator. From there we use the C function itoa() to convert the temp variable to a string and add it to our string.
From here we shove that into either a nearpoints() function or a pcfind() function. In my case I’m running with pcfind() since it seemed to perform a tad bit faster with this. But either method works.
int temp = @cell;string group = "\@cell!=" + itoa(temp);
int find[] = pcfind(0, group, "P", @P, 100, 1);
@cell2 = point(0, "cell", find[0]);
WARNING: Depending on the amount of RAM we have I would suggest not trying this with more than 250k points. This thing will without a doubt crash not just Houdini but your whole computer, given the opportunity. This shit is dangerous.
DANGEROUS LIKE SOME OTTER BACK SCRATCHES
Alright team, here’s where thing get a bit more involved with nodes outside of wrangles. I wish I could wrangle everything, but sometimes it’s better to just relent and let the Houdini gods take the wheel.
For each voronoi cell,
Create a list of neighbor cells from the stored attributes
If there is a list of neighbors from a previous iteration
Append new list to the end of the previous list, while checking for duplicates
Alas it is time to drop into the second fabled For Each Block. Let’s set the end block to run over points, and set the piece attribute to our @cell attribute. The question for us now is, how do we take our cell2 attribute (our neighbor cell), and create a list from all the independent attributes stored on their respective points? The obvious and most straightforward would be to use a loop in VEX to grab each point’s neighbor number, and append it to a detail attribute. Luckily I never go with the straightforward approach, so let’s see if we can’t abuse the point wrangle to give us the result we want without looping.
My first thought was to throw down an Attribute Promote SOP, which is the easiest way to shift data from point level to any other class level, however it doesn’t give use the option to create a list from an attribute.
Let’s do some clever VEX instead:
int nb[];
nb[0] = i@cell2;setdetailattrib(0, "cell2", nb, "append");
This is method is not intuitive, so I understand any confusion that comes from it. In order to use the “append” mode of the functions setdetailattrib, setpointattrib, etc. we are required to first store the data we wish to append, in an array. My theory is that the “append” function matches the array type first, and if it can’t it just defaults to the “set” mode instead. If we don’t first store the data in an array, your detail attribute would likely refer to point(0, “cell2”, npoints(0)); or in layman’s speak, the final point in your VEX loop. Regardless, this will give us an array of all your neighboring cells.
int nb[];
nb[0] = i@cell2;setdetailattrib(0, "cell2", nb, "append");setdetailattrib(0, "cell", @cell, "set");if(len(i[]@nearp) > 0){
setdetailattrib(0, “cell2”, i[]@nearp, "append");
}
The other two things we’ll need to do, is save our current cell to a detail attribute, as well as retain our previous neighbor array, if we have one. Above is the code for that.
Append new list to the end of the previous list, while checking for duplicates
So this is definitely the most standard of array operations, and there are plenty of ways of doing this. I’m not going to go in depth on what my code does, as it is hopefully straight forward:
int temp[] = i[]@cell2;
resize(@cell2, 0);for(int i = 0; i < len(temp); i++){
float find = find(@cell2, temp[i]);
if(sign(find) == -1){
append(@cell2, temp[i]);
}
}
Basically test if the value from the temp array has been added to the cell2 array already, if it hasn’t shove that shit in. Make sure the wrangle is set to details mode as well, otherwise say good bye to your RAM.
The last and final step is to transfer the data back to the points we sampled from initially, and then FEEDBACK (yeezy, yeezy, what’s good).
Before I go any further, I have to say we need to give our vPoints IDs before we transfer the data onto them, so before the main feedback loop, we can use whatever preferred method we have of giving each point an ID.
Now let’s drop down a final wrangle, pipe in our vPoints into Input0 and our A points that we’ve been manipulating into Input1 and slay this dragon.
int nearp[] = detail(1, "cell2");
int cell = detail(1, "cell");
int cell2 = detail(1, "cell2");if(cell != @id){
removepoint(0, @ptnum);
}else{
i[]@nearp = nearp;
@cell = cell;
@cell2 = cell2;
}
This is hopefully the most straightforward chunk of code so far, basically if the current ID isn’t the cell from before, delete the point. If it is the cell we’ve been transforming, add the attributes we’ve been playing with, onto it.
From there we pipe that into the End block of the For Each Block, and then pipe that into the end block of our Feedback Loop. One last thing to note, is that in my test file i’m using the feedback iterations to jitter my surface samples, the A points. That jitter is what allows for us to find new neighbors (or not if we’ve found em all), on each new iteration.
Like last time, I’m not going to cover my methods of meshing this, that portion is already built for you in the file. It’s basically going to connect your current point to its neighbors naively. And then from there we’re going to go through and re-wind those polygons to contain a third neighbor, giving us a full triangle. Neato burrito. However this meshing will only work with the attribute names we’ve specified, if you’re trying to do your own version you will need to adjust the variables accordingly.
Now if we fuse those points together to create full mesh, we can finally compute the dual on it! Here’s what that ultimately looks like.
So fuckin cool.
To recap we literally went from a low res point cloud, to a full triangle mesh, with pretty minimal effort using Voronoi diagrams. We can then compute new discrete Voronoi Diagrams over that mesh which are much easier to work with than methods using magic colors.
That’s pretty incredible. What makes it more incredible is that it doesn’t rely on any surface attributes to be successful, meaning it can be used for tons of applications. Another useful feature is that we can define the total number of points we want out in our final mesh, making it pretty awesome for remeshing with realtime in mind.
Currently the issues with the method are edge loss, meshes tend to minimize with this method, so I’m currently looking into methods of preventing that. Another one is that overall, it’s pretty damn slow compared to Triangulate2D or the Remesh SOP.
I was inspired to build this system while reading through a paper on Laplace Beltrami operators over point clouds. I still have no idea what Laplace Beltrami operators totally are, but it got me to this end result, so I can’t complain. After talking on the Houdini Discord (come join us!) with a few artists about this method, the lovely Petz chimed in and sent me a related algorithm known as the Power Crust Algorithm. Maybe one day I’ll take a stab at that as well.
At this point we’re eleven pages in on Google Docs single spaced so I’ll keep this brief, I hope this was digestible. This blog post is definitely targeted at other technical artists, so I understand how it might be a bit intense if you’re just getting into Houdini. Okay super intense. Just keep at it, and I promise things get easier and easier. The other thing is that I will be posting up a video on Vimeo soon, describing other uses for this, as well as other ways of meshing 3d diagrams using volumes! And like always, hire my ass. Ya boi gotta pay rent. Maybe I’ll start a Patreon :thinking_emoji:
*Just to cover my bases, I have paid for both Houdini and Cinema at this point, and do not advocate piracy in any way (please no sue, thanks).