This is the first installment of a my little tutorial series on mesh to sparse voxel conversion. Conversions between different representations are well described in articles and could be considered trivial. But there are always some challenges involved when trying to implement these kind of things. And it is often the case that reference implementations require a lot of dependencies. So i have tried to compile something that could proves useful into a small an detached code chunk.

Voxelization gone wrong

The first step will be on the topic of converting a triangle mesh into a solid voxel model.

Next we will have a look at how to implement fast marching that produces high quality distance fields. Distance fields require a lot of memory to be accurate, so to the third tutorial will be on the subject of creating adaptive distance fields that are super cool. In the end I will show some code from other guys in our team, where we go full circle and produce triangle and tetrahedral meshes based on the distance fields. So it is quite some work, but it is has proven extremely useful and it is fun to code. And if you are like me, then your also need to see the actual code in order to comprehend all the important details.

In the following i will assume that you have a bit of XP in computer graphics. So the first thing is to load the mesh, and to create an acceleration structure we use to decide where the voxels are. To do this i decided to use a good ol Octree, but with some additional information store. If you are unsure what an octree is you can read more here https://en.wikipedia.org/wiki/Octree.

I guess there is a fair chance that you have used an Octree before, for culling or ray tracing. And guess what! we are going to do exactly that in this tutorial.

Bunny hand subdivided, but unharmed!

Octree voxelization pseudo code.

1. Input your favorite mesh, flattened into a vertex list ( (v0,v1,v2),  (v0,v1,v2).. ), and the voxelsize. Which will determine the number of subdivisions of the Octree.
2. Then calculate the enclosing power of two bounding box. Because we always subdivide each axis into two equal sized parts.
3. Each Octree node is subdivided into 8 children blocks, like seen in the bunny 2d example above.
4. For each of there children we check if there are any intersection between the bounding box and any mesh triangle. This is done using some Thomas Akenine-Möller code 
5. If there are intersection in step 4.) then we either subdivide this INTERNAL node further or if the bounding box is voxel sized, we have a LEAF node.
6. If not intersection we flag the node as EMPTY_LEAF.

So again pretty basic stuff. But to avoid to much triangle aabb intersection checking, we create a triangle index list with all the intersecting triangles that we submit to the recursive subdivision.

To make the octree more useful later on each LEAF is updated with the distance to the minimum distance to the mesh surface. This will prove very handy when we going to create a signed distance field. Because the voxels contains subvoxel distance values.

float fMinDist = numeric_limits<float>::max();
for(uint v = 0; v < kNodeTriangleIndices.size(); v++)
{
  RFVector3f kTriangleVerts[3];
  uint uiTriIndex  = kNodeTriangleIndices[v];
  kTriangleVerts[0] = kMeshVertices[uiTriIndex];
  kTriangleVerts[1] = kMeshVertices[uiTriIndex+1]; 
  kTriangleVerts[2] = kMeshVertices[uiTriIndex+2];
  RFVector3f kNormal = VectorCross(kTriangleVerts[2]-kTriangleVerts[0],kTriangleVerts[1]-kTriangleVerts[0]);
  RFVector3f kClosestPoint = ClosestPointOnTriangle(kTriangleVerts,pkLeafNode->m_kWorldAabb.GetCenter());
  RFVector3f kDelta  = kClosestPoint - pkLeafNode->m_kWorldAabb.GetCenter();
  float fSign = (VectorDot(kDelta,kNormal) < 0.0f)? -1.0f : 1.0f;
  float fDist = VectorLength(kDelta);
  if(fDist < abs(fMinDist)) fMinDist = fDist*fSign;
}
pkLeafNode->m_fDistanceToSurface = fMinDist;

Above is the code that  finds the min dist between intersected triangles and the LEAF bounding box.

So to explain what is going on is that the function ClosestPointOnTriangle(Triangle, pos)[3], returns the closest point inside the triangle. And using this information we create the Delta vector where using a dot product can determine on which side of the triangle the LEAF’s center is. Again stuff needed for the signed distance fun..

This concludes the subdivision, not so hard right?

But to make a solid voxel model we now need to know i the if the EMPTY_LEAF’s are inside or outside the mesh. My solution is as usual to use ray tracing

So basically for each EMPTY_LEAF we shoot a bunch of rays in all directions and determine a voxel coverage value, which tell you how much an EMPTY_LEAF voxel is occlude. And given a use specified Coverage value, we flag a EMPTY_LEAF as inside if the calculated occlusion is above the input coverage value. This is quite clever because the EMPTY_LEAF can consists of many voxels, so we can decide for all of them in one go.

So what remains to be described is the ray tracing of an Octree, which fortunately is extremely simple to do.

bool RayOctreeInterSection(const RFRay3f &kRay, const RFVoxelOctreeNode *pkNode)
{
   if(pkNode == nullptr) return false;
   if(pkNode->m_eNodeType == NT_LEAF)
   {
     return RayBoxIntersection(kRay, pkNode->m_kWorldAabb);
   }
   if(pkNode->m_eNodeType == NT_INTERNAL)
   {
      if(RayBoxIntersection(kRay,pkNode->m_kWorldAabb))
      {
         for(uint i = 0; i < 8 ; i++)
         {
           if(RayOctreeInterSection(kRay,pkNode->m_pkChildren[i]))
              return true;
         }
      }
      else
      {
        return false;
      }
     }
return false;
}

Hey that is it! just some recursive traversal and a bunch of ray box intersections and then we know if something has been hit.

To conclude the Coverage test, we use a Hammersley  uniform spherical sample generator [6] . In the example code 200 rays are shot for each EMPTY_LEAF node to determine coverage.

uint uiCoverage  = 0;
uint uiMisses    = 0;
for(uint j = 0; j < m_kSampleDirections.size(); j++)
{
  RFRay3f kRay(pkNode->m_kWorldAabb.GetCenter(),m_kSampleDirections[j]);
  if(RayOctreeInterSection(kRay,pkOctreeRoot))
  {
    uiCoverage++;
  }
   else
  {
    uiMisses++;
  }
}

float fCoverage = static_cast<float>(uiCoverage)/static_cast<float>(m_kSampleDirections.size());
if(fCoverage >= fMinSolidCoverage)
{
   pkNode->m_fDistanceToSurface = -1.0f; // Inside == Negative
}
else
{
    pkNode->m_fDistanceToSurface = 1.0f; // Outside == Positive
}

Below is the result of an Armadillo voxelization. The green voxels are the leafs and the red and yellow are the empty internals. Where the red are large than voxelsize.

The remaining step is to subdivide the red voxels into smaller and out put the complete voxel data in your favorite format.

That is it for now, below is link to a zipped folder containing the necessary code, please feel free to comment on the article and the code.

Roger over for now Peter!

Click here to download the zipped code.

(Download link updated 12/10 2019) Please contact me if it does not work!

To build the code use cmake. There are no dependencies beside our small RFMath lib which is included. Otherwise it is plain c++ code. I decided not to do a visualization which would blur the code. So this is left to the reader

The demo code loads the armadillo and creates a voxelization based on the above mentioned parameters. It outputs an ply  file with the voxelized mesh as colored vertices. But you could of course easily change this to output to some voxel format of your choice. But i have tried to make it really simple, and to view the result you can download the cool and very handy http://meshlab.sourceforge.net/ 

The sparse voxelization inside meshlab, the black to white are the leaf nodes with the distance values. And blue are voxel sized Empty leafs. The green are the sparse big Empty leaf voxels.

References:

1. http://www.realtimerendering.com/intersections.html  (Mother of all intersection collections)
2. http://cs.lth.se/tomas_akenine-moller (Father of a lot of intersection test )
3. https://en.wikipedia.org/wiki/Octree
4. http://www.gamedev.net/topic/552906-closest-point-on-triangle/ (Link to the closest point code origin)
5. http://meshlab.sourceforge.net/
6. http://holger.dammertz.org/stuff/notes_HammersleyOnHemisphere.html