3D Bitmaps, Minecraft Cubes, and Mesh Winding Numbers
/As a follow-up to my Signed Distance Fields tutorial, a reader asked about how to voxelize a mesh. The MeshSignedDistanceGrid that we computed is a kind of voxelization, but this reader was (I assume) asking about binary voxels - ie the blocky, minecraft-y kind - which can be represented with a Bitmap3 object, where each (i,j,k) entry is true if it is inside the mesh and false if it is outside.
There are several ways to create this Bitmap3 voxelization of a mesh. If you start with the MeshSignedDistanceGrid from the SDF Tutorial, then you can convert it to a binary bitmap like so:
// create SDF
MeshSignedDistanceGrid levelSet = ...
Bitmap3 bmp = new Bitmap3(levelSet.Dimensions);
foreach(Vector3i idx in bmp.Indices()) {
float f = levelSet[idx.x, idx.y, idx.z];
bmp.Set(idx, (f < 0) ? true : false);
}This block creates a Bitmap3 with the same dimensions as the SDF, and then for each index, sets the voxel as "inside" (true) if the distance is negative.
Creating a Minecraft-Style Surface Mesh
If you would like to see what the binary voxelization looks like, you can use the VoxelSurfaceGenerator class to create a Minecraft-style mesh of the voxel faces:
VoxelSurfaceGenerator voxGen = new VoxelSurfaceGenerator();
voxGen.Voxels = bmp;
voxGen.ColorSourceF = (idx) => {
return new Colorf((float)idx.x, (float)idx.y, (float)idx.z) * (1.0f / numcells);
};
voxGen.Generate();
DMesh3 voxMesh = voxGen.Meshes[0];
Util.WriteDebugMesh(voxMesh, "your\\path\\mesh_file.obj");Click to enlarge
The ColorSourceF function I am setting here is used to provide a solid color for each block. This is helpful for visualization, but might also be useful if for example you had some spatial coloring function you wanted to visualize. You can also assign colors to the mesh afterwards, of course.
The last line writes out the mesh to an OBJ file, you'll have to provide your own path. The result is shown to the right, for our standard bunny mesh. This image is a screen shot from Autodesk Meshmixer, with boundary edges shown in blue. Note that the mesh generated by VoxelSurfaceGenerator is actually a bunch of small squares, that are simply adjacent. So, many mesh processing techniques will not work on this mesh.
We could try to weld these borders together, but if we aren't careful this will result in non-manifold topology where more than two faces meet along an edge or at a vertex. We'll cover that in a future tutorial.
Voxelization with a Point-Containment Queries
Generating an SDF is one way to voxelize a mesh, and perhaps the fastest, because we only have to resolve inside/outside near the surface. A fast-sweeping algorithm is then used to fill the rest of space. However, we have some other options, which do have some benefits. Here is an alternative that uses DMeshAABBTree3.IsInside() to set the inside/outside value at each voxel:
DMesh3 mesh = (load mesh...);
DMeshAABBTree3 spatial = new DMeshAABBTree3(mesh, autoBuild: true);
AxisAlignedBox3d bounds = mesh.CachedBounds;
int numcells = 32;
double cellsize = bounds.MaxDim / numcells;
ShiftGridIndexer3 indexer = new ShiftGridIndexer3(bounds.Min, cellsize);
Bitmap3 bmp = new Bitmap3(new Vector3i(numcells,numcells,numcells));
foreach (Vector3i idx in bmp.Indices()) {
Vector3d v = indexer.FromGrid(idx);
bmp.Set(idx, spatial.IsInside(v));
}Basically the same code as above, the main difference is that we have to create a ShiftGridIndexer3 to map between grid coordinates and 3D space. In this example we only need to go from the grid (i,j,k) to 3D (x,y,z), using FromGrid(). However if in your application you want to map from 3D coordinates into the grid, you can use ToGrid(). There are also more advanced Indexer options, for example if your object has a full 3D position (origin + rotation), you can use a FrameGridIndexer3 to do the mapping.
This code uses a smaller grid resolution - 32 - so the result is blockier, as you can see on the right. There are not a lot of benefits to using the mesh IsInside() in this context, however you might find the code above useful if you have some kind of geometry that is hard to get into the level set format.
Voxelization with the Mesh Winding Number
There is a third function we can use to determine inside/outside of a mesh at a point, that is pretty awesome. Here is the code:
spatial.WindingNumber(Vector3d.Zero); // seed cache outside of parallel eval
Bitmap3 bmp = new Bitmap3(new Vector3i(numcells+2, numcells+2, numcells+2));
gParallel.ForEach(bmp.Indices(), (idx) => {
Vector3d v = indexer.FromGrid(idx);
bmp.SafeSet(idx, spatial.WindingNumber(v) > 0.5);
});Ok, the main thing is that we are using DMeshAABBTree3.WindingNumber(), which computes a hierarchical evaluation of the Mesh Winding Number. This is something that was just recently invented, and published in a SIGGRAPH paper in 2013 by Alec Jacobson, Ladislav Kavan, and Olga Sorkine-Hornung. If you have ever used the Polygon Winding Number to check if a point is inside a polygon, this is the same thing, but for 3D meshes. You just compute a simple function over the mesh triangles, and if the mesh is closed, you get an integer if you are inside the mesh, and 0 if you are outside. Kind of amazing.
The function DMesh3.WindingNumber() will do this computation, and for small meshes or single queries (eg like checking if a single 3D point is inside a mesh in a VR user-interface), this is sufficient. However if the mesh is large and/or you are doing lots of queries (like to fill a 64^3 voxel grid), it is too slow. However there is a neat trick to do a hierarchical computation based on a bounding-box tree, which is implemented in DMesh3AABBTree, and that's what I'm using above (see the paper for more details).
Since each voxel is independent, we can trivially parallelize this evaluation. The WindingNumber() function in DMesh3AABBTree does a precomputation the first time it is called, so the first line does one evaluation to seed this cache. Then we use gParallel.ForEeach, which does a multi-threaded iteration over the IEnumerable we provide (here the indices of the grid), and calls the lambda function for each index. Note that we also use Bitmap3.SafeSet() (instead of Set), which internally uses a SpinLock to make sure we aren't writing to the Bitmap3 from multiple threads. This is necessary because internally a BitArray is used, which stores the bits packed into 32-bit integers. Since we cannot read-and-flip bits independently in an atomic operation, we might end up with race conditions without the lock.
Ok, enough of that - why bother with this much more expensive computation? Because it is magic. Below is another example, where I have cut a few large holes in our venerable bunny. Also, that sphere stuck on it's back, is in fact only half a sphere, that is just overlapping - it is not connected to the bunny mesh. Look what happens when we use IsInside() (middle) or the SDF version (right). In both cases, those open boundaries are a disaster.
Now, the interesting thing about the Mesh Winding Number (MWN) is that unlike the binary IsInside(), or the distance-field-based SDF, it is a real-valued computation that is well-defined over space. If the mesh is closed the MWN is an integer, but when the mesh contains holes, the winding number smoothly diffuses around the open boundaries. As a result, we can use a non-integer threshold to determine the range of MWN values we consider to be "inside". In the code above I used > 0.5, which produces a great result shown below-left. This value is important though - if I use 0.8, I get the result on the right - not a catastrophic failure, but there are still clearly big chunks missing.
The Mesh Winding Number is a very powerful tool, and is also extremely simple to implement - my brute-force implementation is about 20 lines of code (the hierarchical version is harder, but not /that/ hard). The hierarchical evaluation is critical if you want to fill a grid. On my 6-core/12-thread machine, in Release build, computing the 64^3 voxelizations above takes a second or two. With the non-hierarchical version, it took well over a minute.
So, go forth and voxelize!
[Update] A reader asked about splitting up a mesh into 64k chunks, which is a hard constraint in the Unity Mesh class. There is a way to do this for any DMesh3, but that's a topic for a future tutorial. The VoxelSurfaceGenerator has a field MaxMeshElementCount, if you set this to 65535, then once the growing mesh hits that many vertices or triangles, a new mesh will be started. The resulting meshes will be accumulated in the VoxelSurfaceGenerator.Meshes List.
[Second Update] If you want to try using the Mesh Winding Number to repair a mesh, my new tool Cotangent uses this algorithm to repair meshes in the Solidify Tool. However it doesn't produce chunky voxels, it produces a smoother surface.