Brute-Force Voxelization

Brute-Force Voxelization

In game dev, nailing performance is like painting a masterpiece – you tweak and refine, but eventually, you have to step back. It’s never quite ‘done,’ merely released. Contrast that with the ‘10% Rule’ of my university professor: if it won’t net you a double-digit boost, why bother? Just run it on a bigger cluster. To most game programmers this sounds like a blasphemy. But even in game development, there are instances where this brute-force approach actually hits the mark. In my case it was voxelization.

Reliable Tool

To give a bit of background, my profesor was working in the field of computational fluid dynamics, so he was more concerned with the accuracy of the physical model rather than performance. In the end it doesn’t really matter if the calculations took 10 or 11 hours if the results are wrong. The 10% rule was not encouraging sloppy code, it was about prioritizing results.

Right now I would benefit from similar result-focused approach. In my destruction tech demo, there are multiple instances where I need to determine how much two irregular shapes are overlapping.

The overlap data will be then used by the physics engine to determine the strength of the bond between 2 objects. There will be thousands of such pairs, so manual setup is impractical. These calculations will be done in editor as a part of destruction authoring tool, and performed on user request, not every frame. This imposes a set of requirements:

  • Ease of implementation: The aim of the tool is to reduce the amount of work, not trade time spent on manual setup for time spent coding.
  • Robustness: It should handle all edge cases without requiring manual correction.
  • Performance: It may take several minutes to complete, I’m not planning to use it in real-time; it will solely function within the editor.
  • Accuracy: Deviations of up to 30% in overlap volume are acceptable.

In essence, the algorithm should mimic human judgment in determining if two meshes overlap sufficiently. While not requiring absolute accuracy, it must avoid false positives – visibly separated meshes should not be classified as overlapping.

Options

I didn’t want to spend days on researching the subject, so I only considered 3 most obvious approaches:

  • AABB overlap
  • Boolean operation on meshes
  • Boolean operations on voxelized meshes

Axis-Aligned Bounding Boxes Overlap

Using Axis-Aligned Bounding Boxes (AABBs) to determine overlap volume is straightforward: it involves finding the AABB for each mesh and calculating the overlap volume. However, this approach can lead to false positives in certain scenarios.

Boolean Operation on Meshes

Boolean operations on meshes are simple in theory. Find the overlapping part and calculate its volume. However, in practice, these operations are known for their unreliability due to issues such as geometric complexity, numerical precision errors, and topological problems.

Boolean Operation on Voxelized Objects

Boolean operations on the voxel objects – voxelize both meshes using the same grid, and then use boolean AND to get the common part. Easy to implement, and robust, provided that the voxelization was performed correctly. Moreover, with sufficiently high resolution, it can yield good accuracy.

Voxelization

While boolean operations on voxelized meshes are reliable, the accuracy of the results hinges on the quality of voxelization. In this particular case the voxelization is binary – per each voxel, or volume pixel, a boolean value is stored – true if the center is inside mesh, false if outside:

Both objects are voxelized using the same grid. If values of a voxel is true for both objects, then the objects overlap in that voxel:

One crucial aspect yet to be clarified is the definition of “inside” and “outside.”

Inside vs Outside

While it may seem immediately apparent whether a point resides inside a shape, translating this concept into an algorithm can be tricky. ‘Even–odd rule’ is helpful here, provided that the shape is closed. This rule determines whether a point lies inside a closed shape by counting the number of times a ray from the point intersects with the shapes’ outline; if the count is odd, the point is inside, and if it’s even, the point is outside.

If the total number of intersections is odd – 1, 3, 5, 7… the point is inside:

If the number is even – 0, 2, 4, 6… the point is laying outside:

Real-life situations are rarely as straightforward as textbook examples. Even when most meshes are designed to be watertight, there are often irregularities like odd shapes, gaps, or overlaps. These irregularities can cause thin, elongated lines to stick out from the voxelized object:

3D Voxelization

The ‘Utah Teapot’ will be the test case, as is a great example of a shape that looks solid from the outside but is tricky to voxelize properly:

Voxelization is an ideal candidate for parallel computing – tens of thousands of simple operations on the same set of data. I will use GPU to speed up the calculations.

First I have to convert the mesh into a form, that will be easy to process. Typically, meshes are stored as lists of vertex positions followed by sets of indices forming triangles. To simplify the code, I’ll embed each vertex’s position directly within its corresponding triangle. All mesh triangles are processed in this manner and then stored in a buffer:

struct Triangle
{
    float3 vertexA;
    float3 vertexB;
    float3 vertexC;
};

StructuredBuffer<Triangle> TriangleBuffer;

With the mesh prepared, the next step is to develop a version of the ‘even-odd rule’ algorithm that operates with an array of triangles. This algorithm will iterate through all mesh triangles for each voxel, checking for intersections:

bool IsInside(float3 origin, float3 direction)
{
    int intersections = 0;
    for (int i = 0; i < TriangleCount; i++)
    {
        if (IntersectTriangle(origin, direction, TriangleBuffer[i]))
        {
            intersections ++;
        }
    }
    return intersections % 2 != 0;
}

For each voxel, a ray is cast from its origin in the direction specified by the user, which remains consistent for all voxels. Ray-triangle intersection checks are carried out using the Möller–Trumbore method:

bool IntersectTriangle(float3 origin, float3 direction, Triangle tri)
{
    float epsilon = 0.000001f;
    float3 edgeAB = tri.vertexB - tri.vertexA;
    float3 edgeAC = tri.vertexC - tri.vertexA;
    float3 pVector = cross(direction, edgeAC);
    float determinant = dot(edgeAB, pVector);
    if (determinant > -epsilon && determinant < epsilon)
    {
        return false;
    }
    float inverseDeterminant = 1.0f / determinant;
    float3 tVector = origin - tri.vertexA;
    float u = dot(tVector, pVector) * inverseDeterminant;
    if (u < 0.0f || u > 1.0f)
    {
        return false;
    }
    float3 qVector = cross(tVector, edgeAB);
    float v = dot(direction, qVector) * inverseDeterminant;
    if (v < 0.0f || v + u > 1.0f)
    {
        return false;
    }
    float t = dot(edgeAC, qVector) * inverseDeterminant;
    if (t > epsilon)
    {
        return true;
    }
    return false;
}

Voxelization Failures

As expected, the results were not without their problems:

I tested several ray directions, and each voxelization exhibited numerous artifacts. To compound the issue, the protruding deformations can intersect with nearby objects. This could potentially result in false-positive overlap checks, precisely the situation I’m trying to prevent.

A quick inspection of the teapot reveals the root of the problem. The mesh is plagued with topological issues such as holes, gaps, and overlaps.

Brute Force Voxelization

When designing the Space Shuttle Orbiter, NASA incorporated five General Purpose Computers (GPC) into its avionics. The idea behind using multiple computers in the Space Shuttle’s flight control system was to ensure redundancy and reliability. By employing several identical computers running in parallel, any failing hardware could be isolated without disruption, as the majority-voting system allowed the operational computers to override the malfunctioning one. The system provided fail-operational, fail-safe capability of the spacecraft’s avionics.

I plan to apply a similar logic to reduce errors in the voxelization algorithm. Each voxel will undergo voting across different directions to decide if it’s inside or outside. Using a ‘majority voting system’, where each direction’s vote counts, should give us a mostly accurate result, with fewer errors.

While using the space shuttle analogy to explain multisampling may seem like a bit of a stretch, it does a great job of illustrating the concept. With such setup I should be able to improve the accuracy by increasing the number of raycasting directions.

Voxelization Directions

Currently, the raycasting directions have been defined by the user. However, as I plan to use an arbitrary number of directions, I’ll need an algorithm to generate them automatically. While this is straightforward in a 2D scenario – simply dividing 360° by the number of directions – it’s less obvious in 3 dimensions.

Luckily, evenly distributing points on a sphere has been tackled before and is known as the ‘Thomson problem’. It describes how electrons would arrange themselves when placed on a unit sphere.

Electrons, being negatively charged, repel each other, leading to an increase in distance from neighboring particles until they reach a point where they can no longer move further. At this stage, we can assume that they are evenly spaced.

Thomson problem iterations

The same procedure can be applied to distribute directions evenly. Here, the end of each direction vector can be likened to an electron that repels other directions. The function takes an array of directions and executes one step of distribution, returning the degree of change in the arrangement. The first direction in an array is treated as stationary.

public static float Distribute(Vector3[] directions, float stepSize, float strength)
{
    if (directions.Length <= 1)
    {
        return 0;
    }
    float maxDelta = float.MinValue;
    for (int i = 1; i < directions.Length; i++)
    {
        Vector3 correction = Vector3.zero;
        for (int j = 0; j < directions.Length; j++)
        {
            if (i != j)
            {
                Vector3 difference = directions[i] - directions[j];
                Vector3 direction = difference.normalized;
                float distance = difference.magnitude;
                float falloff = Mathf.Pow(Mathf.Clamp01((2 - distance) / 2), 2);
                float magnitude = falloff * strength;
                magnitude = Mathf.Min(magnitude, stepSize);
                correction += direction * magnitude / directions.Length;
            }
        }
        Vector3 oldDirection = directions[i];
        directions[i] += correction;
        directions[i] = Vector3.Normalize(directions[i]);
        maxDelta = Mathf.Max(maxDelta, Vector3.Distance(oldDirection, directions[i]));
    }
    return maxDelta;
}

With this function we can start with any random arrangement of directions and iterate the Distribute step until the directions are evenly spaced.

List<Vector3> GenerateDistributed(int count)
{
    List<Vector3> directions;
    // Generate random directions, start with one fixed fixed direction
    directions = GenerateRandom(count, Vector3.up);

    // Distribute the directions until convergence or 1000 iterations
    int iteration = 0;
    float maxDelta = 0.0f;
    while (iteration < 1000 && maxDelta >= 0.001f)
    {
        maxDelta = Distribute(directions, 0.01f, 100000);
        iteration++;
    }
    return directions;
}

Here, the algorithm is presented one step per frame, making it easier to follow along with the process:

It works well for any number of directions, although the convergence may be slow in some cases.

Result

Although sampling multiple directions is a crude approach to solving the problem, it is undoubtedly effective. Increasing the number of samples can largely mitigate voxelization errors. Remarkably, satisfactory results for the teapot test case can be achieved with just 5 directions.

Saved Time

At every step of the process, I opted for the path of least resistance—from checking every voxel in the domain to pushing away the directions. I chose the simplest and probably most computationally intensive algorithms available.

Making the connection graph for the stone wall asset takes about a second. Although it’s not perfect, this step is only done once per asset. With tens of assets expected, a faster algorithm could potentially save a few minutes at best. However, developing such an algorithm would likely require about a week. Not exactly a good trade.

Conclusions

In cases when you are not trying to fit within 16.6 milliseconds, but within the deadline, simplicity often trumps sophistication. Completing the job “the right way” isn’t always more valuable than completing it on time.

I could have invested weeks or even months in extensive research and meticulous code polishing, only to potentially end up with a system that nobody would appreciate or use. Instead, I opted for shortcuts and hacks, and surprisingly, they did the job. The voxelization system is functional, was quick to implement, and now I have the time left to spend on other problems.