Merging Meshes with Signed Distance Fields
/
Left: lots of overlapping spheres, YOUR SLICER's NIGHTMARE
Right: a single continuous surface OF ULTIMATE PRINTABILITY
[Update July 6, 2018] My new tool Cotangent exposes the mesh-to-SDF-to-mesh operation I describe below in the Solidify Tool, if you want to try it without writing C# code! [/EndUpdate]
In this tutorial I'll show you how to use a few of the tools in geometry3Sharp to create solid approximations to an input mesh. We'll create a Signed Distance Field (SDF) approximation for an input DMesh3, and then use MarchingCubes to extract a new mesh from the SDF. These basic techniques can be used for an enormous number of cool things, from voxelizing and mesh repair (similar to the Make Solid tool in Autodesk Meshmixer), to data structures for spatial queries, and even for 3D modeling things like Mesh Booleans.
As a teaser, by the end of this tutorial, I'll have explained how to turn that mess of overlapping spheres on the right into a closed surface that you could safely send to any 3D printing software.
By the way, if you would like to experiment with g3Sharp, it's now available as an official NuGet package - you can find it by searching for "geometry3Sharp" in your favorite package manager.
I'll start with a basic code sample and then we'll step through it:
DMesh3 mesh = StandardMeshReader.ReadMesh("c:\\demo\\bunny_solid.obj");
int num_cells = 128;
double cell_size = mesh.CachedBounds.MaxDim / num_cells;
MeshSignedDistanceGrid sdf = new MeshSignedDistanceGrid(mesh, cell_size);
sdf.Compute();
var iso = new DenseGridTrilinearImplicit(sdf.Grid, sdf.GridOrigin, sdf.CellSize);
MarchingCubes c = new MarchingCubes();
c.Implicit = iso;
c.Bounds = mesh.CachedBounds;
c.CubeSize = c.Bounds.MaxDim / 128;
c.Bounds.Expand(3 * c.CubeSize);
c.Generate();
DMesh3 outputMesh = c.Mesh;
StandardMeshWriter.WriteMesh("c:\\demo\\output_mesh.obj", c.Mesh, WriteOptions.Defaults);That's it. If you run this code on this solid bunny mesh, then open the input and output meshes, you'll see that they look quite similar. The SDF version is a bit smoother in some places, and it has more triangles so the shading is different. But if you overlay the two, you'll see that they solid is a very close approximation (right image).
Lets step through the code. After loading in the mesh, I first decide on num_cells. This defines the density of the grid we will compute the SDF on. Larger numbers mean better shape approximation, but also more memory and computation time. Internally, the SDF will be based on a dense 3D grid of floats - think of a solid block of Minecraft voxels. Each block is a cube cell_size wide.
Next we create the MeshSignedDistanceGrid object for the mesh and compute the result, which comes out in sdf.Grid. Then we create a DenseGridTrilinearImplicit based on this grid. This class will use trilinear interpolation to turn the discrete grid values (ie the blocky cubes / voxels) into a continuously-varying 3D scalar field. So, based on our Signed Distance Grid, we have created a Signed Distance Field.
We can call iso.Value(Vector3d) at any point in space, and it will return the (approximate) signed distance to the (approximated) mesh surface. If the distance is negative, then that point is inside the surface. Positive is outside, and if the value is zero, then the point is on the surface (this surface is known as the iso-surface). Of course we rarely get to exactly zero, but if we had point A just inside the surface (ie negative), and point B just outside (positive), then we know that at some point on this line, the function will evaluate to zero. So, we can do root-finding along this line, using something like Bisection or Newton's method, to converge on the zero value (wikipedia).
This is exactly how we will get a mesh back out of our SDF. Remember, at this point our SDF is completely decoupled from the original mesh. All we have is the function iso.Value(). So what we are going to do is fill a bounding-box with smaller boxes, and evaluate the SDF at the corners of all the boxes. Then when we find a box where some corners are inside and some are outside, we know the surface cuts through that box. We'll do root-finding along the box edges to (approximately) find the zeros, and then make a patch of triangles. This is the famous Marching Cubes algorithm. Paul Bourke has an excellent page with more information, and you can find a PDF of the original paper here - compared to a modern SIGGRAPH paper it is incredibly readable.
In g3Sharp the MarchingCubes class implements this method, and you can give it an arbitrary function to surface via the Implicit member. This way of representing a 3D shape - as a scalar function over space - is often called an Implicit Surface, hence the naming. We also have to provide a bounding box to search inside of (ie to fill with smaller cubes), and the CubeSize we wish to use. Note the small expansion of the bounds - if you don't do this, the most extreme edges of the shape might be clipped off (3 cells is probably overkill).
After calling Generate(), the mesh is built and we can write it out. Easy!
You're not impressed? Ok, how about this one. Here's the same bunny with an intersecting sphere stuck on it, and I did a bit of sculpting to create a self-intersection (download). If you wanted to turn this into a solid non-self-intersecting shell (say, to 3D print it), well, you're in trouble. These are hard problems to solve. But not for an SDF - run this file through the code above, and the output mesh is a closed, manifold shell. The trade-off for this simplicity is that we have to accept some resampling artifacts around the sharp edges.
(click to enlarge)
In the examples above, what is "inside" vs "outside" is defined by something called the winding number. We'll explore this concept more in a future post. But, basically, in this context it means that for any given point in space, we can count up how many "times" the point is inside the input surface. So, points inside both the sphere and the bunny, have winding number 2. Similarly, points inside that bit that is self-intersecting (which you can't see) also have winding number 2 - they are inside the bunny "twice". Points outside have winding number 0.
But what about cavities on the inside of the surface? Well, the winding number depends on the orientation of the mesh. If we flip a mesh inside-out, then the points inside it have negative winding numbers. In the SDF mesher, we will define any point with positive winding number as inside. This means we can use inside-out meshes to define holes, and also to do Boolean subtraction. Here's an example below. I am using Meshmixer to view the output meshes, and in Meshmixer the red/white stripey pattern means you are looking at the "back-side" of the mesh surface. I cut away in the SDF+MC mesh to show that there is also a fully-enclosed hole.
click to embiggen
In the description above, I mentioned two parameters - the cell_size we passed into MeshSignedDistanceGrid, and the MarchingCubes.CubeSize. The cell_size of the Signed Distance Grid defines how accurate our approximation of the input mesh can be. Even if we mesh at a very high resolution, if the grid doesn't capture the shape well, we'll have approximation error.
In the image below, I set num_cells to 16, 32, 64, 128, and 256, and used a CubeSize of 128. At the lower resolutions we clearly don't have enough values in the grid to capture the shape. However it is quite interesting to mesh a low-resolution grid at high resolution. The trilinear interpolation produces a smoothly-varying patch inside each cell, but you can clearly see the discontinuities at the cell boundaries, where the interpolation is only C0 continuous. At the higher 128 and 256 resolutions, the grid is as or more accurate than the Marching Cubes triangulation, so you don't see any difference.
click for big
Varying the CubeSize is not so nice at lower resolutions - again I went from 16-256 below, this time with num_cells=256 for the SDF grid. The only reason you might use lower-resolution marching cubes is to reduce memory (Marching Cubes will produce huge meshes very quickly!). At the highest resolution below, you can start to see the original triangles - clearly our SDF approximation is quite accurate in this case! However even at this resolution, the sharp edge around the border of our intersection is not particularly clean. This is caused by both the SDF and the Marching Cubes, but even if it were a nice, smooth sharp edge in the SDF, Marching Cubes will not be able to capture it.
bigbig
So, to capture shape details we want to use high MC resolutions, but then we end up with huge meshes. What to do? Use another tool in our mesh processing toolbox, of course. The Reducer class, which I described in a previous tutorial, can help us get back to reasonable mesh resolutions. In this case you only need a few lines:
Reducer r = new Reducer(outputMesh);
r.ReduceToTriangleCount(50000);In the 256-SDF-grid, 256-MC-grid case above, the initial MC mesh has about 450k triangles. The code above will reduce this mesh to 50k triangles. This takes about 3 seconds on my 4-core desktop computer. This is actually more time that it takes to generate the 450k mesh in the first place! The result, shown in the middle on the right, is much lower-density but clearly there is still a lot of redundant geometry. If we go further, down to 10k triangles (far right), the mesh starts to get better.
You'll also notice that in the 10k-triangles version, the sharp edge around the intersection has started to get a bit cleaner. I have (experimentally) found that using a slight modification to the Reducer setup will do an even better job at recovering these sharp edges. Instead of reducing to a specific triangle count, the Reducer also supports reducing to a target edge length. The following code:
r.ReduceToEdgeLength(2*c.CubeSize)resulted in the mesh on the right, which has a very crisp sharp edge around where the sphere was subtracted. This doesn't always work, but it does work sometimes!
So, we can combine a multiple overlapping meshes into a single SDF solid, mesh it at a crazy high resolution, and then use Mesh Simplification to get it back to a closed shell we could actually use in other software. And maybe even get sharp edges out. What to do? Generate some shapes, of course! Geometry3Sharp has lots of different kinds of Mesh Generators built in - Spheres, Cylinders, Boxes, Tubes, Surfaces of Revolution, and more. Look in the mesh_generators/ folder to find all the generator classes. You can do procedural shape generation in two steps - first generate a bunch of small closed meshes, then combine them all and send through the code above. It literally is that easy.
As an example, here is a bit of code that generates a new mesh from an input mesh, by adding a sphere for each vertex, and a rectangular box along each edge:
Sphere3Generator_NormalizedCube gen = new Sphere3Generator_NormalizedCube() { Radius = sphere_radius, EdgeVertices = 5 };
DMesh3 sphereMesh = gen.Generate().MakeDMesh();
DMesh3 latticeMesh = new DMesh3();
MeshEditor editor = new MeshEditor(latticeMesh);
foreach ( int vid in mesh.VertexIndices() ) {
DMesh3 copy = new DMesh3(sphereMesh);
MeshTransforms.Translate(copy, mesh.GetVertex(vid));
editor.AppendMesh(copy);
}
foreach ( Index4i edge_info in mesh.Edges() ) {
Vector3d a = mesh.GetVertex(edge_info.a), b = mesh.GetVertex(edge_info.b);
Frame3f f = new Frame3f((a + b) * 0.5, (b - a).Normalized);
editor.AppendBox(f, new Vector3f(box_width, box_width, (b - a).Length*0.5));
}I set the sphere_radius and box_diam to something appropriate for the scale of my mesh, and ran it on a bunny reduced to 512 triangles, with the SDF-grid cell_size and Marching Cubes CubeSize both set to 512. This crazy resolution is required to capture the fine details. After about 30 seconds, the set of overlapping meshes on the left is turned into the single solid on the right:
You really want to click on this one!
Hopefully this gives you some ideas =)
One caveat about Marching Cubes should be mentioned, and is illustrated in the low-res MC image further up the page. I mentioned that we figure out the mesh patch that should be inside each "cube" based on whether the corners of the cube are inside or outside the surface. Most of the time, this works well, but there are cases where it is ambiguous. Currently, the Marching Cubes implementation in g3Sharp does not properly handle these cases. It is mostly a problem when the shape varies much more rapidly than the MC mesh. In our case these failures will leave small holes in the mesh, as certain triangles will be discarded (in other implementations these triangles would produce non-manifold geometry, but DMesh3 doesn't allow that). We'll be fixing this in the future.
To give credit where credit is due, the super-fast Mesh to SDF conversion I have implemented in g3Sharp is based on a C++ implementation by Christopher Batty, which is on github. And he tells me this code was a cleaned-up version of an initial implementation by Robert Bridson. My C# version has been extensively refactored and modified so that some steps can be multi-threaded. I also added parity-counting to properly handle overlapping shells.