algorithmic modeling for Rhino

For another project I have been working with metaball scalar fields, and I thought it would be fun to have a go at meshing them using Paul Bourke's marching tetrahedra methods.  I know that Daniel Piker has done a little of this before more generally, and that there is an excellent implementation of the marching cubes algorithm as a part of Millipede...this is really just more specifically set up for generating metaball meshes.

On its own, the marching tetrahedra algorithm as outlined by Paul Bourke produces some odd results for metaballs...basically, the points that drive the mesh are assigned along the edges of the tetrahedra according to a linear interpretation of the iso value...what this means is that if you have an "ideal" iso value say of 0.29, and the sample points at each end of your tetrahedron edge return values of 0.2 and 0.3, the point will be located along that edge at 10% of the length away from the point with the 0.3 sample.  In general, this works fine, except (from what it appears to me at least...could be wrong here) that the metaball field itself isn't really linear, and the tetrahedra used for sampling are carved from cubes such that the lengths of each edge differ that there are a few spikes and dips here and there.  Sampling resolution has a lot to do with these results as well.  I address this issue here by embedding laplacian smoothing as a part of the process, and also by an additional field value sampling run for each mesh vertex after the initial mesh has been constructed.

So the inputs are your metaball centers (as a list), their radii (also as a list, which is the location at which their isosurface value contribution is 0), the iso level threshold, the grid resolution for sampling (here a lower number produces a higher resolution...basically is the rough size of the grid in the x, y and z directions), and the number of smoothing passes run after the mesh is generated.  The outputs include the metaball meshes as well as the bounding box that's used, as well as the centers of the sample cells (both of these last are superfluous but perhaps interesting).  I'd suggest starting with a low resolution (higher res value) for sketching, as it works pretty quickly...then if you want, you can increase the the sampling density (and incur the resulting exponential slowdown).

I've used a couple of tricks for reducing the number of sampled points (nothing so sophisticated as an octree subdivision), but there are also definitely inefficiencies in points shared between sampling cells having their fields calculated multiple times.  There's really plenty of room for improvement...this is just a rough pass really, but I figured some people might have some fun with it and it seems that every now and again the issue of making metaball meshes comes up, and I figured that it'd be useful to have a quick way of building (relatively) clean meshes laying about.

Also, probably lots of bugs everywhere...feel free to point them out...and the code's pretty ugly, but poke around, take it, do what you will with it.

Views: 6200

Tags: mesh, metaballs, vb


You need to be a member of Grasshopper to add comments!

Join Grasshopper

Comment by David Stasiuk on February 21, 2015 at 6:34am

Hi Rodion- sorry I missed your post on Wednesday! I would really stay away from's super heavy and totally unnecessary, as the marching tetrahedrons algorithm organizes all of your points perfectly for creating your own mesh faces to begin with. You really don't need to build surfaces at all, or get points...

I don't use Python, so I can't describe the syntax there, but the C# syntax for creating a mesh from those points would be something like this:

First, create your final mesh object:

Mesh MyFullMesh = new Mesh();

Then, run your marching tetrahedrons algorithm, and for each Triangle output (which should give you three points) do this:

Mesh NewTriangle = new Mesh();






MyFullMesh.Vertices.CombineIdentical(true, true);

Sorry that it's late...hope it helps! 

Comment by Rodion Kalistsuk on February 21, 2015 at 6:06am

Hello Again,

So I managed to do the 1st 2 points. Apparently, in Weaverbird there is a component Laplacian smoothening. So here are my results:

This is fine for me already, I guess. The only thing, is that the meshing command (done in the 2nd result) is done via Rhino plugin: MeshFromPoints. I don't know yet how to access it in grasshopper or in Iron Python, so that I would have a script, where I don't need anything to do manually...

Thank You David for Your help !

Comment by Rodion Kalistsuk on February 18, 2015 at 4:22pm

Hey David,

I didn't have the time to work on that code. Now I am going to do it within these couple of days.

I wanted to clarify what I have to do in those points which You suggested:

1) build up a full mesh rather that bunch of surfaces, and use Rhinocommon to combine identical vertices, and rebuild the vertex normals

As the output of the code I have are the points that are generated by the triangular planes, I can just use these points with a plugin called MeshFromPoints, which generates a mesh from the points you give it. However, the problem of bad mesh is still there.

As You may also observe, yes there are probably points that are duplicated, however, there are points that are really close to each other, that produce the irregularity of the mesh.

Now from what You have written, I can delete the duplicates, OK, but what do You mean by rebuild the vertex normals? Could You please explain it in detail?

2) run a couple rounds of laplacian smoothing on the mesh to better distribute your vertices (for each vertex, make it equal in location to the average of its neighbours)

Here, I also have a question. Ok, I can find the average of the neighboring points by math, but how do I determine within scripting the neighboring points ?

3) create a line normal to each vertex roughly the length of your sampling grid and test the endpoints of it against your scalar field formula, and then do one final linear interpolation between those two points for your vertex

This one is again not clear what to do. Could You please explain ?

Comment by David Stasiuk on February 6, 2015 at 1:37am

Hi Rodion-

I think you're confusing marching cubes with marching're right that marching cubes has 256 cases, but the lookup you're using for marching tetrahedrons is correct. There are just the 8 cases, so you're actually doing it right here (scroll down on this page, and you'll see a separate subset all about marching tetrahedrons The benefit to using marching tetrahedrons is exactly this: that the number of possible "cuts" through the tetrahedron are dramatically smaller in number than those through a cube.

However, I have found that also what you're seeing that the linear interpolation creates some odd distortions (which is why I went ahead and later did the marching cubes implementation). Some of this comes from the density of the sampling grid: the more dense, the fewer distortions.

What I would suggest, if you want a (relatively) quick way to improve this outcome:

1) build up a full mesh rather that bunch of surfaces, and use Rhinocommon to combine identical vertices, and rebuild the vertex normals

2) run a couple rounds of laplacian smoothing on the mesh to better distribute your vertices (for each vertex, make it equal in location to the average of its neighbours)

3) create a line normal to each vertex roughly the length of your sampling grid and test the endpoints of it against your scalar field formula, and then do one final linear interpolation between those two points for your vertex.

This should give you a smoother mesh for sure.

But good work getting this far! 

Comment by Rodion Kalistsuk on February 5, 2015 at 6:03pm

Dear David,

With the help of my friend, who is a much more experienced programmer than I, after some tweaking and playing around with the output and seeing what it does, we managed to receive the Metaballs. So yes, the output of that code are points, generated by the Marching cubes algorithm.

this is the picture:I didn't apply any meshing. I just created the separate surfaces between the sets of 3 points defining a triangle. Just to visualize the result in a better way.

As You can see, the pattern of points in certain places is rather distorted. Could You please help with solving this problem?

In the code, which is discussed earlier, there is a comment and some lines:

def PolygoniseTri(g, iso, v0, v1, v2, v3):
triangles = []

# Determine which of the 16 cases we have given which vertices
# are above or below the isosurface

triindex = 0;
if g.val[v0] < iso: triindex |= 1
if g.val[v1] < iso: triindex |= 2
if g.val[v2] < iso: triindex |= 4
if g.val[v3] < iso: triindex |= 8

return trianglefs[triindex](g, iso, v0, v1, v2, v3)

But according to Paul Bourke's article, there should be 256 cases of triangular intersections, so:

cubeindex = 0;    
if (grid.val[0] < isolevel) cubeindex |= 1;
if (grid.val[1] < isolevel) cubeindex |= 2;
if (grid.val[2] < isolevel) cubeindex |= 4;
if (grid.val[3] < isolevel) cubeindex |= 8;
if (grid.val[4] < isolevel) cubeindex |= 16;
if (grid.val[5] < isolevel) cubeindex |= 32;
if (grid.val[6] < isolevel) cubeindex |= 64;
if (grid.val[7] < isolevel) cubeindex |= 128;

Is this the reason, or it is because of something else?

This is rather important, because it is an important aesthetic aspect of our Final Project.
Thank You in advance

Comment by David Stasiuk on February 4, 2015 at 4:43am

Sounds like you're on the right track, Rodion. I really haven't dug into this particular code very closely, but from what it appears to me, yes, the PolygoniseTri should give you the vertex locations for any new mesh faces derived from that particular sampling tetrahedron. And it looks to me like it's correctly cycling through each tetrahedron in each cube:

  1. triangles.extend(PolygoniseTri(grid,isolevel,0,2,3,7))  
  2.     triangles.extend(PolygoniseTri(grid,isolevel,0,2,6,7))  
  3.     triangles.extend(PolygoniseTri(grid,isolevel,0,4,6,7))  
  4.     triangles.extend(PolygoniseTri(grid,isolevel,0,6,1,2))  
  5.     triangles.extend(PolygoniseTri(grid,isolevel,0,6,1,4))  
  6.     triangles.extend(PolygoniseTri(grid,isolevel,5,6,1,4)) 

So in grasshopper, this is where you would define a new mesh face, using the points indicated from this function (looks like they return either one or two triangles depending on the intersection state of the tetrahedron, which is right). By iterating through your entire sample grid, you'll create your continuous mesh.

I have found that the easiest way to do it is to create a new Mesh for each sample tetrahedron, and then to append them to a complete mesh as you go along. Hope this'll be cool to make a Python implementation for GH!

Comment by Rodion Kalistsuk on February 3, 2015 at 4:26pm

Dear David,

I understood the most of the code, but I have some questions, if You could help me figuring them out, it would be great!

I understand, that the program starts with the function: def readdata(f=lambda x,y,z:x*x+y*y+z*z,size=5.0,steps=11):
where this part:f=lambda x,y,z: is an anonymous function, with the help of which, we calculate the values at each point. It calculates it with respect to the points generated in the function: def lobes(x,y,z): which is a mathematical function for a certain shape, which is shown in the picture in the web page of the code.

Then, in the main() function, the program flow is continued by assigning the previously calculated values to the vertices of the rectangles, with which the 3D domain is discretized. So we receive something similar to the picture in the code web page, under the chapter 'Solution'.

After this, each rectangle is generated as an object 'grid' of the class 'Gridcell', assigning the previously calculated values and vertices as 2 arrays to each 'grid' object. 

Now, in the function PolygoniseTri(grid,isolevel,0,2,3,7the program determines, which points are below the isosurface according to the algorithm that is described in . But, it is not clear, why is the function called only with 4 vertices, so then they are checked only with 4 conditions: 
  1.  if g.val[v0] < iso: triindex |= 1  
  2.  if g.val[v1] < iso: triindex |= 2  
  3.  if g.val[v2] < iso: triindex |= 4  
  4.  if g.val[v3] < iso: triindex |= 8
Meaning that only vertices 0 2 3 7 are checked, for this example. What about the other 4 ?
The return statement of this function is: return trianglefs[triindex](g, iso, v0, v1, v2, v3) On line 186, while trianglefs = {7:t0708,8:t0708,9:t0906,6:t0906,10:t0A05,5:t0A05,11:t0B04,4:t0B04,12:t0C03,3:t0C03,13:t0D02,2:t0D02,14:t0E01,1:t0E01,0:t000F,15:t000F} we see it is a 16 optioned result.
Ok, moving on. As far as I understand, return trianglefs[triindex](g, iso, v0, v1, v2, v3) returns the necessary call of the one of the functions like def t0708(g, iso, v0, v1, v2, v3): which correctly linearly interpolates the position of the iso surface on an edge. So basically it creates the triangles. Or more likely its points.

So as I understand the output are points of the resulting Metaball surface? Im not sure about this.

If it is possible, give Your comments on this.

Comment by Rodion Kalistsuk on February 2, 2015 at 3:01pm

Thank You very much David,

I will look into this implementation and try to adapt it!

Comment by David Stasiuk on February 2, 2015 at 1:44am

This looks like it is a pretty reasonable Python implementation, based on the same source that I used for mine:

You'd have to adapt it to Rhinocommon, of course, but all of the key ingredients are there.

Comment by Rodion Kalistsuk on February 1, 2015 at 6:36am

Hello David,

Yea. One of the conditions of the course final project is to use Python, since we were learning scripting the geometry in Python.

I guess we are going to try to go around the problem, but in future I guess I will look into these implementations, just for myself. 



Search Grasshopper


  • Add Photos
  • View All

© 2015   Created by Scott Davidson.   Powered by

Badges  |  Report an Issue  |  Terms of Service