# GPU Mesh Voxelizer Part 5: Filling the inner voxels

In this article, we’ll fill the inside voxels of our voxelized mesh. Doing this will complete our GPU-based mesh voxelizer. To do this, we’ll implement the Möller-Trumbore algorithm for ray/triangle intersection. If you haven’t seen the previous parts, you can read them here: part1part2part3, and part4.

## How to fill the inside voxels

I’ll describe the algorithm before we begin. For every voxel, we cast a ray in any direction to find the first triangle intersection. If the dot product of the ray and the triangle normal is greater than zero, it implies that we’re exiting the mesh. That means we must’ve started inside the mesh. I’m not sure if this algorithm covers every possible edge case. However, from what I can tell, it seems to work well enough. Initially, I used a different method, but I had to swap it out at the last minute when I discovered voxels that were mistaken to be inside the mesh.

Now that we understand the goal, how can we achieve it? Well, we can’t perform raycasts on the GPU yet, so let’s start there. Let’s implement the Möller-Trumbore algorithm.

## The Möller-Trumbore algorithm

What’s the Möller-Trumbore algorithm? It’s a fast way to compute ray/triangle intersection, invented by Tomas Möller and Ben Trumbore. I won’t explain all the math behind it, but I’ll leave some resources at the end of this article if you’re interested. Since I’m not going to explain the math, I’ll go ahead and drop the function here.

``````// VoxelizeMesh.compute
bool IntersectsRayTriangle(float3 rayOrigin, float3 rayDirection, Triangle tri, out float hitDistance, out float3 hitNormal)
{
const float epsilon = 0.0000001;

float3 ab = tri.b - tri.a;
float3 ac = tri.c - tri.a;

float3 h = cross(rayDirection, ac);
float a = dot(ab, h);

if (a > -epsilon && a < epsilon) return false;

float f = 1.0 / a;
float3 s = rayOrigin - tri.a;
float u = f * dot(s, h);
if (u < 0.0 || u > 1.0) return false;

float3 q = cross(s, ab);
float v = f * dot(rayDirection, q);
if (v < 0.0 || u + v > 1.0) return false;

float t = f * dot(ac, q);
float3 bc = tri.c - tri.b;
hitNormal = cross(ab, bc);
hitDistance = t;
return t > epsilon;
}``````

To use this function, we’ll write a new kernel in the compute shader. That’s because I decided to break up the process into two steps. First, we calculate the shell of the mesh, and then we fill the volume. Breaking it up made it easier to write and debug, but it would be interesting to combine them and test the performance difference. So let’s get to it.

## The Fill Volume kernel

Let’s add the new kernel to the `VoxelizeMesh.compute` file. By the way, if you don’t want to follow along, I’ll link the project at the end of this article.

``````#pragma kernel VoxelizeMesh
#pragma kernel FillVolume
...
{
if (id.x >= _GridWidth || id.y >= _GridHeight || id.z >= _GridDepth)
return;

const int voxelIndex = id.x + _GridWidth * (id.y + _GridHeight * id.z);
if (_Voxels[voxelIndex].isSolid > 0.0)
return;
}
...``````

Here’s the boilerplate out of the way. We added a new kernel, and in that kernel, we exit early if we’re outside the bounds of the voxel grid or if the voxel is already solid. So now we have to start raycasting against triangles. Typically we would do a broad phase first to rule out as many triangles as possible quickly, then raycast against what’s left. However, for simplicity, I’m raycasting against every triangle for each voxel. So that means there’s plenty of room for optimization. Luckily GPUs are very fast, which means we can get away with this for now. So, let’s iterate over every triangle and find the closest hit.

``````...
const float cellSize = _CellHalfSize * 2.0;
float3 rayOrigin = mul(_WorldToLocalMatrix, float4(
id.x * cellSize + _CellHalfSize + _BoundsMin.x,
id.y * cellSize + _CellHalfSize + _BoundsMin.y,
id.z * cellSize + _CellHalfSize + _BoundsMin.z, 1.0));

float3 rayDirection = float3(0, 0, 1);
float minDistance = 10000000.0;
float3 minDistanceNormal;
bool didIntersect = false;

for (int i = 0; i < _TriangleCount; i += 3)
{
Triangle tri;
tri.a = _MeshVertices[_MeshTriangleIndices[i]];
tri.b = _MeshVertices[_MeshTriangleIndices[i + 1]];
tri.c = _MeshVertices[_MeshTriangleIndices[i + 2]];
float3 outNormal;
float distance;
if (IntersectsRayTriangle(rayOrigin, rayDirection, tri, distance, outNormal))
{
if (distance < minDistance)
{
minDistance = distance;
minDistanceNormal = outNormal;
didIntersect = true;
}
}
}
...``````

We start by calculating the origin of the ray. Remember, we’re running one thread per voxel. With one ray per thread, each ray is at the center of a voxel. Then, we continue by choosing a ray direction and initializing our minimum distance. By the way, it doesn’t matter what direction you pick because when you’re inside the mesh, every direction will lead out before going back in. Then we check every triangle for a ray intersection and store the normal of the closest hit. Then let’s finish by checking the dot product of the normal and the ray. If the dot product is greater than 0, then our ray points in more-or-less the same direction as the triangle normal, which implies it’s exiting the mesh. If the ray is leaving the mesh, then we can assume this voxel is inside the mesh.

``````...
if (didIntersect && dot(minDistanceNormal, rayDirection) > 0.0)
_Voxels[voxelIndex].isSolid = 1.0;
}``````

Here’s the entire function, so it’s easier to read.

``````[numthreads(4,4,4)]
{
if (id.x >= _GridWidth || id.y >= _GridHeight || id.z >= _GridDepth)
return;

const int voxelIndex = id.x + _GridWidth * (id.y + _GridHeight * id.z);
if (_Voxels[voxelIndex].isSolid > 0.0)
return;

const float cellSize = _CellHalfSize * 2.0;
float3 rayOrigin = mul(_WorldToLocalMatrix,
float4(id.x * cellSize + _CellHalfSize + _BoundsMin.x,
id.y * cellSize + _CellHalfSize + _BoundsMin.y,
id.z * cellSize + _CellHalfSize + _BoundsMin.z, 1.0));
float3 rayDirection = float3(0, 0, 1);
float minDistance = 10000000.0;
float3 minDistanceNormal;
bool didIntersect = false;

for (int i = 0; i < _TriangleCount; i += 3)
{
Triangle tri;
tri.a = _MeshVertices[_MeshTriangleIndices[i]];
tri.b = _MeshVertices[_MeshTriangleIndices[i + 1]];
tri.c = _MeshVertices[_MeshTriangleIndices[i + 2]];
float3 outNormal;
float distance;
if (IntersectsRayTriangle(rayOrigin, rayDirection, tri, distance, outNormal))
{
if (distance < minDistance)
{
minDistance = distance;
minDistanceNormal = outNormal;
didIntersect = true;
}
}
}

if (didIntersect && dot(minDistanceNormal, rayDirection) > 0.0)
_Voxels[voxelIndex].isSolid = 1.0;
}``````

That wraps up the GPU side of things. All we have to do is call our kernel from our C# script. Over in `VoxelizedMesh.cs` add this block to the end of `VoxelizeMeshWithGPU()`, after dispatching the first kernel.

``````var volumeKernel = _voxelizeComputeShader.FindKernel("FillVolume");
_meshVerticesBuffer);
_meshTrianglesBuffer);