Merging Meshes with Signed Distance Fields

Left: lots of overlapping spheres, YOUR SLICER's NIGHTMARE Right: a single continuous surface OF ULTIMATE PRINTABILITY

Left: lots of overlapping spheres, YOUR SLICER's NIGHTMARE
Right: a single continuous surface OF ULTIMATE PRINTABILITY

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);

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);

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.


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);

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:


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));
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. 

Mesh Simplification with g3Sharp

Recently a user posted a github issue asking for a mesh simplification example. It just so happened that I had recently finished writing an implementation of Garland and Heckbert's Quadric Error Metric (QEM) Simplification algorithm. If you want to learn more about this technique, the original papers and several later articles are available on Michael Garland's Website, and are very readable. I will give the broad strokes below, but first here is an example of what we are talking about - automatic reduction of a bunny mesh from 13,000 to 500 triangles:


An easy way to Simplify or Reduce a mesh (which is the terminology I use, because...well I don't have a good reason, but it's how I think of it!) is to iteratively collapse edges in the mesh. Each collapse removes one edges and the two (or one at the boundary) triangles connected to that edge, as well as one vertex. Just repeat until you hit your target triangle count, and you're done. Easy! 

Except, which edge should you collapse first? And, when you collapse the edge, should you just keep the existing vertex position? In most cases it would at least make more sense to move the remaining vertex to the edge midpoint. But could we do better? This is where QEM comes in. Basically, the Quadric Error gives us a way to (1) "score" edges, so we know which edge collapses will have the smallest impact on the shape, and (2) predict a position for the new vertex that will minimize the score after the collapse. 

Ultimately, the Quadric Error is a measurement of the sum-of-distances from a point to a set of planes. If you think of the ring of triangles around a vertex, each triangle makes a plane.  The distance from the vertex to that plane is zero. But if we start moving the vertex, this distance increases. So we can "score" an a vertex movement (like an edge collapse) by measuring all the distances from the new vertex to all the input planes (of the original triangles). Each point-plane distance measurement can be expressed as a matrix multiply with the point, and since Ap+Bp = (A+B)p, we can combine all the error measurements into a single matrix multiplication!

Even better, after an edge collapse, we can think of that vertex as having a total error measured relative to the input planes of both vertices (still just one matrix). So as we do sequential collapses, we accumulate all the plane-distance-functions of the original mesh triangles that were (at some point in the past) connected to the remaining vertex. And it's still just one matrix. In the '99 paper linked from the site above [PDF], Garland showed how the QEM error is in some sense a measure of surface curvature, and produces "optimal" triangulations, under a reasonable definition of optimal. Amazing!

But, ok, you're probably just here for the code, right? To use the g3sharp implementation you need to get your mesh to be a DMesh3, see my previous tutorial for how to do that. Then you just need a couple lines to create a Reducer object and run the simplification:

DMesh3 mesh = load_my_mesh_somehow();
Reducer r = new Reducer(mesh);

In the code above, 500 is the target triangle count. This takes a fraction of a second to run. I have not done extensive profiling, but on a relatively fast machine I can reduce a 350k mesh to 100k in 2 to 3 seconds. So, it's pretty fast.

Mesh Validity Checking

Most of the mesh processing algorithms in g3Sharp require that the input meshes be manifold. This term has lots of meanings, in this context it means that at minimum each mesh edge is connected 1 or 2 triangles. Inside DMesh3, edge triangles are stored in an Index2i object, so DMesh3 can't even represent a non-manifold edge. In most cases we also require that there are no bowtie vertices, which are vertices connected to disjoint sets of triangles. The simplest example of this would be two triangles connected only at one vertex (hence the name "bowtie"). 

If these conditions are not met, then some mesh operations - like an edge collapse - can produce undefined results, result in exceptions, and so on. So, to test that your mesh is internally consistent, you can use the DMesh3.CheckValidity() function. You can configure this function to have different behavior if the mesh is not valid, ie it can assert, throw an exception, just return false, and so on. By default it will consider bowtie vertices to be invalid, but you can also configure this behavior with the bAllowNonManifoldVertices argument. 

Note that this function also does very thorough testing of the internal mesh data structures to make sure everything is consistent. So, it is relatively expensive, and you probably do not want to be calling it all the time in production code. One exception is when loading meshes from a file, in that case you really should check the meshes after reading, it can save you a lot of headaches trying to track down garbage-in/garbage-out type problems!

Preserving Boundaries

In some cases, for meshes with boundaries you might need to preserve the boundary loops exactly. For example if you are actually reducing a sub-region of a larger mesh, or a mesh split along UV boundaries. This just takes a few more lines, before your call to ReduceToTriangleCount():

r.SetExternalConstraints(new MeshConstraints());
MeshConstraintUtil.FixAllBoundaryEdges(r.Constraints, mesh);

The MeshConstraints object in the code above is in fact a very powerful facility, that allows you to constrain much more than just boundary edges. But that is a topic for a future tutorial. You can poke around the data structure, or the MeshConstraintUtil helper class, to find out more. The images below compare reducing a bunny-with-boundary to 500 triangles without (left) and with (right) a preserved boundary. 


Here is a closeup of the boundary around the bunny's front paws. You see on the left that there are significantly shorter edges along the boundary loop, because it has been exactly preserved. However, you might also note if you look closely (or click to enlarge) that on the front-left paw that there are some thin sliver triangles. This is a current limitation of boundary preservation - it may result in a bit of ugly stuff at the border. This will hopefully be improved in future updates.


Project to Target

Finally, one last thing that you might want to do when Simplifying a mesh. By default, the "new" vertex position after an edge collapse is computed by minimizing the QEM error for that vertex. Compared to something like edge midpoints, this produces nicer shapes and actually results in the algorithm running much faster. However in some cases you may require that the vertex positions lie on the original mesh surface. This is also supported. First you build a spatial data structure, like the DMeshAABBTree3 we built in the Spatial Query tutorial, and then set that as a "Projection Target" for the Reducer. This causes the vertex positions to be mapped to the nearest points on the input mesh surface. Here is the code:

DMeshAABBTree3 tree = new DMeshAABBTree3(new DMesh3(mesh));
MeshProjectionTarget target = new MeshProjectionTarget(tree.Mesh, tree);
r.ProjectionMode = Reducer.TargetProjectionMode.Inline;

The last line is optional. This causes the Reducer to compute the projection each time it wants to evaluate the QEM error for that vertex. This is "more correct" but many of these vertices will eventually be discarded, so the work is in some sense wasted (projections are expensive). If you leave this line out, then the projection is computed after the reducer is finished, for just the vertices that were ultimately kept. The image below compares no projection (left) and with inline projection (right), overlaid on the original surface (dark grey). The simplified meshes don't actually look very different, but you can see that on the right, most of the mesh is "inside" the original surface, while on the left it is roughly half-and-half inside/outside.


One thing to keep in mind with Projection is that for thin parts, the projection can easily end up on the "wrong" side. So, for most simplification problems you probably don't need it. Which is great because it's quite a bit slower!

Now, go save some triangles! Next up is Remeshing, which works in much the same way, only it's also much more complicated...

Mesh Creation and Spatial Queries with g3Sharp

Welcome to the first post on the gradientspace blog! The first of many, I hope. The posts here will mainly be tutorials on how to use the various gradientspace open-source libraries, and, eventually, interactive 3D tools. 

In this first entry, I will answer a user question, which was filed as Issue # 2 in the geometry3SharpDemos project. The gist is that the user would like to use geometry3Sharp to construct a mesh and do some raycasting against the mesh surface. 

The first problem is how to construct a DMesh3 object from lists of vertex x/y/z coordinates, triangle indices, and in this case also normals (which might not be necessary for the user's problem). This is not very hard but it comes up so often in my own coding that I decided to add a new utility function that makes this construction a one-liner:

DMesh3 mesh = DMesh3Builder.Build(vertices, triangles, normals)

This DMesh3Builder.Build() function is written using C# generics, and internally it does type interrogation to figure out what the input buffers are and cast them to the correct types. So, the vertices and normals arguments could be a float[] array, a List<Vector3f>, or any other generic IEnumerable of <float>,<double>,<Vector3f> or <Vector3d> type. Similarly triangles can be an int[] array or any other IEnumerable<int> or <Index3i>.

This uber-function is not necessarily the most efficient. Internally it basically does this:

    DMesh3 mesh = new DMesh3(MeshComponents.VertexNormals);
    for ( int i = 0; i < NumVertices; ++i )
        mesh.AppendVertex(new NewVertexInfo(vertices[i], normals[i]));
    foreach ( Index3i tri in triangles )

The NewVertexInfo type has additional constructors for other cases, such as vertex colors and UVs. Note that you need to bitwise-or in additional flags (eg MeshComponents.VertexColors) in the constructor, or use the functions like DMesh3.EnableVertexColors(), to allocate these other internal data structures before you can add colors.

After you create a mesh like this, it is a good idea to check that all the internal data structures are consistent. In some cases AppendTriangle() will throw Exceptions if there is a problem, but we do not exhaustively check that the mesh is well-formed on construction because those checks are expensive. Instead, you can call DMesh3.CheckValidity() to do this. This function takes a FailMode argument which determines whether it throws, asserts, or returns false when a problem is found.

(If you do find problems, fixing them might be difficult - I recommend trying Autodesk Meshmixer for now...)

Basic Mesh File I/O

After you have constructed a mesh as above, you might want to see what it looks like. You can do this by exporting the mesh to disk and opening it in a mesh viewer, like the aforementioned Meshmixer. The code to write out a single mesh is a somewhat-convoluted one-liner:

    IOWriteResult result = StandardMeshWriter.WriteFile(path,
            new List<WriteMesh>() { new WriteMesh(mesh) }, WriteOptions.Defaults);

OBJSTL, and OFF formats are supported. For STL the default is ASCII, but if you want a smaller binary STL you can configure this in the WriteOptions data structure, long with many other standard and format-specific options. 

If you would like to read a mesh from disk, you can use the StandardMeshReader class. This currently can read OBJSTL, and OFF formats. It is possible to register additional readers yourself using the MeshFormatReader interface. The simplest way to read a mesh is a one-liner:

DMesh3 mesh = StandardMeshReader.ReadMesh(path)

This works for most cases but if your file contains multiple meshes, or you want to get error feedback, or configure read options, you have to use the more verbose method:

    DMesh3Builder builder = new DMesh3Builder();
    StandardMeshReader reader = new StandardMeshReader() { MeshBuilder = builder };
    IOReadResult result = reader.Read(path, ReadOptions.Defaults);
    if (result.code == IOCode.Ok)
        List<DMesh3> meshes = builder.Meshes;

For OBJ format we can also read the materials, but you have to load the texture images yourself. This is somewhat complicated, perhaps a topic for a future post.

Spatial Data Structure Queries

The next part of the Issue asks how to make a spatial data structure, to do efficient ray-intersection queries. Currently g3Sharp only supports Axis-Aligned Bounding Box (AABB) trees. It just takes two lines to set one up:

DMeshAABBTree3 spatial = new DMeshAABBTree3(mesh);

If the mesh is large this might take a few seconds, but the result is a spatial data structure that has many query functions. For example we can compute a ray-cast like so:

Ray3d ray = new Ray3d(origin, direction);
int hit_tid = spatial.FindNearestHitTriangle(ray);

Of course the ray might miss, so we have to check the resulting triangle ID:

    if (hit_tid != DMesh3.InvalidID) {
        IntrRay3Triangle3 intr = MeshQueries.TriangleIntersection(mesh, hit_tid, ray);
        double hit_dist = origin.Distance(ray.PointAt(intr.RayParameter));

Generally when a query returns a vertex, triangle, or edge index, you should test it against DMesh3.InvalidID to check if the query actually found anything. 

DMeshAABBTree3 also supports nearest-point queries, which are very useful in lots of applications. Here is a the standard code to find the nearest point on a mesh to an input point:

    int near_tid = spatial.FindNearestTriangle(point);
    if (near_tid != DMesh3.InvalidID ) {
        DistPoint3Triangle3 dist = MeshQueries.TriangleDistance(mesh, near_tid, point);
        Vector3d nearest_pt = dist.TriangleClosest;

Those are the two queries I use most often, but there are a few others, like FindAllHitTriangles(Ray3d), which finds all ray/triangle intersections, and TestIntersection(Triangle3d), which tests a triangle for intersection with the mesh.

DMeshAABBTree3 also supports a point inside/outside query, using IsInside(point). This function only works if the mesh is closed. In addition, the current implementation is not the most efficient (it uses FindAllHitTriangles() and then counts crossings).

To check for intersections between two meshes,  you can use TestIntersection(DMeshAABBTree3 otherTree). This is more efficient than testing each triangle separately because it descends the bounding-box hierarchies recursively. This function also take an optional Func<Vector3d, Vector3d> TransformF argument, which allows you to apply a transformation to the second mesh without actually modifying its vertex positions. If your meshes are in the same coordinate system you can just pass null for this argument. However if you are, for example, trying to compute intersections between meshes in Unity that have hierarchical transforms above them, then sending in a suitable transform function that maps one into the space of the other can simplify your code.

Finally, if you would like to implement your own spatial queries that can take advantage of the DMeshAABBTree3 spatial decomposition, you can use the internal TreeTraversal class by replacing the box and triangle test functions, and passing your instance to DoTraversal(). See the code comments for more info about how this works.

So, now you know how to load a DMesh3 or construct one from scratch, create an AABB Tree for it, and compute things like raycasts and nearest-points. This is enough to do some pretty interesting procedural geometry generation, if you think about it...