I wanted to use the GPU-based mesh voxelizer we previously wrote to implement an automatic skinning algorithm. We’ll use a heat diffusion algorithm to radiate heat through the voxels to tie bones to vertices. This article focuses on the implementation of the volumetric heat diffusion algorithm. By the way, I’ll post a link to the code on Github at the end. I’ve used the code from the end of the previous project as a starting point.

Automatic skinning through heat diffusion

This article is inspired by an older post by David Rosen about automatic mesh skinning through volumetric heat diffusion. I’ll link that post later, but I’ll quickly explain the idea. If we radiate heat from a bone in our mesh and distribute it through our voxel representation, we can use that heat as reasonably accurate automatic skinning weights. Since the heat can only pass through solid voxels, we have a more precise distance calculation than if we had used the Euclidean distance, which could travel through empty air.

We’ll start by finding the voxel the occupies the position of the current bone. Then, we’ll go through each voxel and set its heat value to the average of its neighbours. Doing a single pass won’t distribute the heat very much, so we’ll run this kernel multiple times until we’re happy with the result. To accomplish all this, we’ll have to make some changes first. Let’s get started.

Adding heat to the Voxel struct

First thing’s first, we’ll need to store each voxel’s current heat value. To do this, we’ll upgrade the voxel struct on both the shader side and the C# side. At the same time, let’s add a field to remember which voxel is currently emanating heat. We’ll use this information to skip this voxel when distributing the heat, keeping it at the highest possible heat level. So if you have the code from the previous project, open VoxelCommon.cginc and add some fields.

// VoxelCommon.cginc
struct Voxel
{
    float3 position;
    float isSolid;
    float heat;
    float isBone;
};

Next, open up VoxelizedMesh.cs and change the Voxel struct on the C# side. At the same time, we’ll have to modify the stride of the Voxel buffer. If you forgot, the stride tells the system how long a single element in a buffer is. So, since we’re adding two floats, we have to increase the stride by the size of two floats.

// VoxelizedMesh.cs
public struct Voxel
{
    Vector3 position;
    float isSolid;
    float heat;
    float isBone;
}
...
void VoxelizeMeshWithGPU()
{
	...
	
	if (resizeVoxelPointsBuffer || _voxelsBuffer == null || !_voxelsBuffer.IsValid())
	{
    		_voxelsBuffer?.Dispose();
    		_voxelsBuffer = new ComputeBuffer(xGridSize * yGridSize * zGridSize, 6 * sizeof(float));
	}
	
	...
}

So now we have a place to store the voxel heat. Before we continue, we need to zero out the heat and the isBone values. We’ll write a new initialization kernel and run it before proceeding with the rest of the calculations.

Clearing the Voxel heat

Let’s open VoxelizeMesh.compute and add a new kernel. As mentioned, this kernel will clear the heat and isBone values.

#pragma kernel ClearBoneHeat
...
[numthreads(4,4,4)]
void ClearBoneHeat(uint3 id : SV_DispatchThreadID)
{
    if (id.x >= _GridWidth || id.y >= _GridHeight || id.z >= _GridDepth)
        return;

    const int voxelIndex = id.x + _GridWidth * (id.y + _GridHeight * id.z);

    _Voxels[voxelIndex].heat = 0.0;
    _Voxels[voxelIndex].isBone = 0.0;
}

This kernel is straightforward. Each thread will grab a voxel and zero it out. Hopefully, you remember how to dispatch the kernel from C#, but if you need help, you can look at the complete project at the end of this post. The next step is to pick a bone and determine which voxel it’s occupying so we know from where we’ll radiate heat.

Finding the starting bone voxel

We’ll write a new kernel that iterates through every voxel and checks if it contains the currently selected bone. By the way, in the engine, a bone is just a pivot point. Later we’ll see why this is an issue, but right now, let’s focus on finding the starting voxel. So, we can already convert from a position to a voxel AABB, but we can’t check if an AABB contains a point yet. Let’s look at the SetStartBone kernel, and then we’ll fill in what’s missing.

float3 _StartBonePosition;
[numthreads(4,4,4)]
void SetStartBone(uint3 id : SV_DispatchThreadID)
{
    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 < 1.0)
        return;

    AABB voxelAabb;
    voxelAabb.extents = float3(_CellHalfSize, _CellHalfSize, _CellHalfSize);
    voxelAabb.center = _Voxels[voxelIndex].position + voxelAabb.extents;

    if (AabbContainsPoint(voxelAabb, _StartBonePosition))
    {
        _Voxels[voxelIndex].heat = MAX_HEAT;
        _Voxels[voxelIndex].isBone = 1.0;
    }
}

As usual, we’ll run one thread per voxel. I’m assuming that each bone is inside the mesh, so if a voxel isn’t solid, we exit early. Otherwise, set up the voxel’s AABB and check if the bone is inside. If it is, set the heat to MAX_HEAT and isBone to 1. By the way, MAX_HEAT is a value I defined in VoxelCommon.cginc, currently it’s set to 100. Now we’re missing the AabbContainsPoint function. As it turns out, that function is straightforward. All we have to do is check if our point is greater than the minimum position and less than the maximum position of the AABB.

bool AabbContainsPoint(AABB aabb, float3 p)
{
    float3 min = aabb.center - aabb.extents;
    float3 max = aabb.center + aabb.extents;

    return min.x <= p.x && max.x >= p.x && min.y <= p.y && max.y >= p.y && min.z <= p.z && max.z >= p.z;
}

With the compute shader portion out of the way, we can focus on the C# side. Since we did this several times already, I’ll go through it rather quickly. First, replace the MeshFilter _meshFilter field with a SkinnedMeshRenderer. We have to assume that we’re dealing with models that have bones now; whether it’s skinned or not doesn’t matter. Using a SkinnedMeshRenderer gives us access to the mesh’s bones. That means we can iterate through them, setting each one as the starting bone position and calculating the bone heat. In the example project, I pick a bone from the inspector rather than iterating through all the bones. I grab the bone through the SkinnedMeshRenderer bones array.

_voxelizeComputeShader.SetVector("_StartBonePosition", _skinnedMeshRenderer.bones[_boneIndex].position);

Afterwards, I dispatch the kernel like all the others. If you want to see the complete code, check out the project linked at the end. Now let’s distribute the bone heat.

Distributing Bone Heat

To distribute the bone heat, we’ll implement another kernel in VoxelizeMesh.compute. Just like the others, we’ll dispatch one thread per voxel. Each thread will select its voxel and each of its neighbours. Then, as long as the neighbour is solid and within the grid’s bounds, we add its heat to the current voxel’s heat and finally divide by the number of available neighbours we found. In other words, we’ll set each voxel’s heat to the average heat of its neighbours. Here’s the code.

[numthreads(4,4,4)]
void DistributeBoneHeat(uint3 id : SV_DispatchThreadID)
{
    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 < 1.0)
        return;

    if (_Voxels[voxelIndex].isBone > 0.0)
        return;

    const int neighbour1 = (id.x + 1) + _GridWidth * (id.y + _GridHeight * id.z);
    const int neighbour2 = (id.x - 1) + _GridWidth * (id.y + _GridHeight * id.z);
    const int neighbour3 = id.x + _GridWidth * ((id.y + 1) + _GridHeight * id.z);
    const int neighbour4 = id.x + _GridWidth * ((id.y - 1) + _GridHeight * id.z);
    const int neighbour5 = id.x + _GridWidth * (id.y + _GridHeight * (id.z + 1));
    const int neighbour6 = id.x + _GridWidth * (id.y + _GridHeight * (id.z - 1));

    float enableNeighbour1 = _Voxels[neighbour1].isSolid * (id.x + 1 >= _GridWidth ? 0.0 : 1.0);
    float enableNeighbour2 = _Voxels[neighbour2].isSolid * ((int)(id.x - 1) < 0 ? 0.0 : 1.0);
    float enableNeighbour3 = _Voxels[neighbour3].isSolid * (id.y + 1 >= _GridHeight ? 0.0 : 1.0);
    float enableNeighbour4 = _Voxels[neighbour4].isSolid * ((int)(id.y - 1) < 0 ? 0.0 : 1.0);
    float enableNeighbour5 = _Voxels[neighbour5].isSolid * (id.z + 1 >= _GridDepth ? 0.0 : 1.0);
    float enableNeighbour6 = _Voxels[neighbour6].isSolid * ((int)(id.z - 1) < 0 ? 0.0 : 1.0);

    const float heat =
        _Voxels[voxelIndex].heat +
        _Voxels[neighbour1].heat * enableNeighbour1 +
        _Voxels[neighbour2].heat * enableNeighbour2 +
        _Voxels[neighbour3].heat * enableNeighbour3 +
        _Voxels[neighbour4].heat * enableNeighbour4 +
        _Voxels[neighbour5].heat * enableNeighbour5 +
        _Voxels[neighbour6].heat * enableNeighbour6;

    float weightSum =
        enableNeighbour1 +
        enableNeighbour2 +
        enableNeighbour3 +
        enableNeighbour4 +
        enableNeighbour5 +
        enableNeighbour6;

    _Voxels[voxelIndex].heat = clamp(heat / weightSum, 0.0, MAX_HEAT);
}

As I mentioned previously, we’ll run this kernel multiple times per frame. Every time we run it, it’ll spread the heat a little further, so it’s up to you to decide how many passes will provide the best distribution. Here’s a sample of the heat distribution.

• 

• 

• 

Results?

So I have to admit, I haven’t finished the algorithm because I haven’t applied the resulting heat back to the bone weights. However, at this point, I’m skeptical that we’ll see good results in the end. The problem is that bones in the engine have a pivot point but no defined endpoint. In a DCC tool like Blender, we could use the bone’s head and tail to calculate a bounding box. We could use the bounding box to find all the voxels that touch the bone, which would give us a more accurate heat distribution. Instead, we only have a single voxel at one end of the bone emanating heat. So, unfortunately, I don’t believe we’ll be able to generate reliable skinning weights.

So it’s likely that this experiment won’t work in the end. I’ll keep pursuing a solution, but don’t expect to see another post on the matter unless I hit a breakthrough. Despite that, I think we still learned a lot. If you know how we can overcome the flaws, please reach out because I would love to hear about it.

Explore the complete project here on GitHub. If you’d like to read the original article by David Rosen, check it out here. If you enjoy reading these articles, join my mailing list to be notified when the next one is released.