GPU Mesh Voxelization 2026

2026-02-28 21:00 UTC gpt-5-2-thinking Open in ChatGPT ↗

If you mean “correct” = conservative (no cracks/holes) surface voxelization + correct inside/outside (solid) classification on the GPU, the newest mesh-voxelization-specific paper I can find (as of early 2026) is:

Delgado Díez et al., 2024 (Eurographics / EG Digital Library): “Optimizing Surface Voxelization for Triangular Meshes with Equidistant Scanlines and Gap Detection” — a single compute-pass surface voxelizer that traverses triangles with parallel, equidistant scanlines and explicitly detects gaps to avoid missing voxels. citeturn0search28turn0search7

That said, for “provably conservative / watertight in practice” GPU voxelization, most modern systems still build on the same core idea:

Practical “most-correct” GPU recipe in 2026 (industry-standard core)

1) Conservative surface voxelization (guarantee coverage)

Use conservative rasterization (hardware if available, otherwise expand triangle footprint) and pick the dominant axis per triangle to minimize holes.
This is the backbone described in classic “conservative voxelization” work citeturn0search17 and still used in newer real-time systems/papers that cite the Crassin-style pipeline citeturn0search6turn0search11.
A clear, modern-ish GPU write-up of this rasterization-style approach is also summarized in TU Wien’s voxelizer PDF (dominant axis + conservative rasterization) citeturn0search8.

2) Solid voxelization (correct interior)

After you have a conservative shell, you need an inside/outside fill:

For watertight meshes, the classic single-pass GPU solid voxelization approach is still a key reference point (and the basic parity/span idea remains common) citeturn0search14.
In practice today: people often do surface voxelization → per-column parity fill (or flood fill from the boundary in a bounded volume). Newer papers/systems still follow this architecture even when the surface step changes. citeturn0search6turn0search11

When to use which “latest” approach

If you want latest surface-voxelization research detail: look at Delgado Díez et al. 2024 (scanlines + gap detection, compute-pass oriented). citeturn0search28turn0search7
If you want maximum correctness in real-time GPU pipelines: conservative rasterization + robust fill (parity/flood), i.e., the conservative voxelization lineage citeturn0search17turn0search8 plus a solidification step citeturn0search14.
If your input meshes are not watertight / “triangle soup”, you’ll usually need a robust preprocessing/repair or robust voxelizer (often CPU-side, sometimes hybrid). Ultraliser is an example framework emphasizing robust voxelization from imperfect inputs (though it’s not presented as “the newest GPU-only algorithm”). citeturn0search2

One quick check so I can point you to the best match

Are you targeting:

surface voxels only (shell), or
solid voxels (filled interior), and do you need it to handle non-watertight meshes?

If you answer that in one line, I’ll recommend the most appropriate current GPU pipeline (DX12/Vulkan/CUDA) and the best-matching references.

Below is (1) a breakdown of what the paper’s shader is doing, and (2) a Unity/HLSL compute shader that mirrors the paper’s functions (CalcPrimaryAxis, SortVertices, CalcFlatAdvance, PointInEdge, Scanline (Amanatides & Woo variant + bias), GapDetection, InsideTriangle, writeOutput), organized in a way that’s practical in Unity.

Everything here is based directly on the paper’s pseudocode listings and descriptions. citeturn0view0turn1view0turn1view1turn1view2turn1view3turn1view4

1) Paper algorithm breakdown (as code steps)

A. Per-triangle initialization

Primary axis selection = least dominant axis of the triangle, i.e., the axis where abs(normal) is smallest (steepest plane gradient in that axis). This axis is treated as “vertical”. citeturn1view1
Sort vertices by height along primary axis so that v0 is highest, v2 lowest, v1 middle. citeturn1view1turn1view0
Flat advance direction = projection of the primary axis onto the triangle plane, then drop the primary component and keep only the sign of horizontal components. This direction drives the scanline traversal bias & “opposite voxel” choice. citeturn1view1

B. Interior traversal via equidistant scanlines

For each integer “floor” h from dot(v2,pAxis) to dot(v0,pAxis):
- Find the segment of intersection between the triangle and the plane dot(P,pAxis)=h by computing two points on triangle edges (PointInEdge).
- Voxelize that segment with a 3D DDA (Amanatides & Woo style) (Scanline). citeturn1view0turn1view2

C. Gap detection (the paper’s key addition)

Each visited voxel checks exactly one neighbor in the current “floor”: its opposite voxel along flatAdvance.
It tests whether the triangle plane intersects the shared edge between the voxel and its opposite within the current floor (distance t <= 1) and whether that intersection point lies inside the triangle (InsideTriangle). If so, it marks the opposite voxel too. citeturn1view3turn1view4

D. Edge processing

After interior scanlines, voxelize triangle edges with the same Scanline function to catch gaps near scanline endpoints and edge-slope cases. citeturn1view0turn1view4

2) Unity compute shader implementing the paper’s functions

This is a single .compute file. It writes a binary surface voxelization into a RWStructuredBuffer<uint> (0/1 per voxel). Use atomics to avoid races.

Assumptions (match the paper’s intent):

Triangle vertices are already in voxel grid space, where 1 unit = 1 voxel, and integer coords align with voxel corners.

You want surface voxels (binary). Solid fill is a separate pass (not in this paper).

Volume dims are _Dim (x,y,z). Bounds are [0..Dim-1] voxel coordinates.

// Unity HLSL compute port of:
// "Optimizing Surface Voxelization for Triangular Meshes with Equidistant Scanlines and Gap Detection"
// Delgado Díez et al., 2024. (Core functions from Listings 1-4)

#pragma kernel ClearVolume
#pragma kernel VoxelizeTriangles

// -----------------------------
// Data
// -----------------------------

struct Tri
{
    float3 v0;
    float3 v1;
    float3 v2;
};

StructuredBuffer<Tri> _Tris;
uint _TriCount;

RWStructuredBuffer<uint> _Voxels;   // length = _Dim.x * _Dim.y * _Dim.z
int3 _Dim;

// Optional: coordinate transform (if you feed object/world space)
// float4x4 _WorldToVoxel; // not used in core listing port

// -----------------------------
// Utilities
// -----------------------------

uint LinearIndex(int3 p)
{
    return (uint)(p.x + p.y * _Dim.x + p.z * _Dim.x * _Dim.y);
}

bool InBounds(int3 p)
{
    return (p.x >= 0 && p.y >= 0 && p.z >= 0 &&
            p.x < _Dim.x && p.y < _Dim.y && p.z < _Dim.z);
}

// writeOutput(voxel): set binary occupancy
void WriteOutput(int3 voxel)
{
    if (!InBounds(voxel)) return;

    uint idx = LinearIndex(voxel);

    // Set to 1 (binary). InterlockedExchange is enough here.
    // If you store bitfields later, change to InterlockedOr.
    InterlockedExchange(_Voxels[idx], 1);
}

float3 CalcNormal(float3 a, float3 b, float3 c)
{
    return normalize(cross(b - a, c - a));
}

// -----------------------------
// Paper functions
// -----------------------------

// CalcPrimaryAxis: axis with smallest abs(normal) component (least dominant axis).
float3 CalcPrimaryAxis(float3 normal)
{
    float3 an = abs(normal);

    // Pick smallest component
    if (an.x <= an.y && an.x <= an.z) return float3(1, 0, 0);
    if (an.y <= an.x && an.y <= an.z) return float3(0, 1, 0);
    return float3(0, 0, 1);
}

// SortVertices: descending order along primary axis
void SortVertices(inout float3 v0, inout float3 v1, inout float3 v2, float3 pAxis)
{
    float h0 = dot(v0, pAxis);
    float h1 = dot(v1, pAxis);
    float h2 = dot(v2, pAxis);

    // We want v0 highest, v2 lowest.
    // Simple sorting network for 3 values.
    if (h0 < h1) { float3 t=v0; v0=v1; v1=t; float th=h0; h0=h1; h1=th; }
    if (h1 < h2) { float3 t=v1; v1=v2; v2=t; float th=h1; h1=h2; h2=th; }
    if (h0 < h1) { float3 t=v0; v0=v1; v1=t; /* heights not needed further */ }
}

// Listing 2: flatAdvance = sign( proj * (1 - pAxis) )
// proj = pAxis - normal * dot(normal, pAxis)
int3 CalcFlatAdvance(float3 pAxis, float3 normal)
{
    float3 proj = pAxis - normal * dot(normal, pAxis);

    // Remove vertical component (the component where pAxis==1)
    float3 horizMask = (1.0.xxx - pAxis); // (0,1,1) etc
    float3 flat = proj * horizMask;

    // Convert to sign in {-1,0,1}
    int3 s;
    s.x = (flat.x > 0.0) ? 1 : ((flat.x < 0.0) ? -1 : 0);
    s.y = (flat.y > 0.0) ? 1 : ((flat.y < 0.0) ? -1 : 0);
    s.z = (flat.z > 0.0) ? 1 : ((flat.z < 0.0) ? -1 : 0);
    return s;
}

// PointInEdge(vA, vB, h): point where dot(pAxis, P) == h
float3 PointInEdge(float3 vA, float3 vB, float3 pAxis, float h)
{
    float a = dot(vA, pAxis);
    float b = dot(vB, pAxis);
    float denom = (b - a);

    // Degenerate edge w.r.t pAxis: return vA (safe-ish)
    if (abs(denom) < 1e-8) return vA;

    float t = (h - a) / denom;
    return lerp(vA, vB, t);
}

// InsideTriangle(point, v0,v1,v2): robust-ish edge test in plane using the triangle normal.
// Accepts points very close to edges (epsilon).
bool InsideTriangle(float3 p, float3 a, float3 b, float3 c, float3 n)
{
    // Ensure consistent orientation with n
    float3 e0 = b - a;
    float3 e1 = c - b;
    float3 e2 = a - c;

    float3 c0 = cross(e0, p - a);
    float3 c1 = cross(e1, p - b);
    float3 c2 = cross(e2, p - c);

    float d0 = dot(c0, n);
    float d1 = dot(c1, n);
    float d2 = dot(c2, n);

    // Allow small negative due to fp
    const float eps = 1e-5;
    return (d0 >= -eps && d1 >= -eps && d2 >= -eps) ||
           (d0 <=  eps && d1 <=  eps && d2 <=  eps); // handle flipped winding
}

// Listing 4: GapDetection
void GapDetection(int3 voxel, int3 flatAdvance, float3 pAxis, float3 normal,
                  float3 v0, float3 v1, float3 v2)
{
    // L0 is the lower corner of the edge connecting voxel to its opposite 12-neighbour.
    // In GLSL listing: L0 = voxel + vec4(greaterThan(flatAdvance, 0), vec4(0.0f));
    float3 L0 = (float3)voxel;
    L0.x += (flatAdvance.x > 0) ? 1.0 : 0.0;
    L0.y += (flatAdvance.y > 0) ? 1.0 : 0.0;
    L0.z += (flatAdvance.z > 0) ? 1.0 : 0.0;

    // P0 is any point on plane, they use top vertex v0
    float3 P0 = v0;

    float denom = dot(pAxis, normal);
    if (abs(denom) < 1e-8) return; // nearly parallel, skip

    // distance along pAxis to intersection
    float distance = dot(P0 - L0, normal) / denom;

    // intersection point
    float3 point = L0 + pAxis * distance;

    // only if intersection occurs within current floor (<= 1 voxel along pAxis)
    if (distance <= 1.0 && distance >= 0.0 && InsideTriangle(point, v0, v1, v2, normal))
    {
        int3 oppositeVoxel = voxel + flatAdvance; // sign(flatAdvance) == flatAdvance here
        WriteOutput(oppositeVoxel);
    }
}

// Listing 3: Scanline (Amanatides & Woo variant + bias opposing flatAdvance)
void Scanline(float3 start, float3 end,
              int3 flatAdvance, float3 pAxis, float3 normal,
              float3 v0, float3 v1, float3 v2)
{
    float3 d = end - start;

    // step = sign(d)
    int3 step;
    step.x = (d.x > 0.0) ? 1 : ((d.x < 0.0) ? -1 : 0);
    step.y = (d.y > 0.0) ? 1 : ((d.y < 0.0) ? -1 : 0);
    step.z = (d.z > 0.0) ? 1 : ((d.z < 0.0) ? -1 : 0);

    float3 fStart = floor(start);
    float3 fEnd   = floor(end);

    // steps = int(dot(floor(end)-floor(start), step))
    int stepsCount = (int)round((fEnd.x - fStart.x) * (float)step.x +
                                (fEnd.y - fStart.y) * (float)step.y +
                                (fEnd.z - fStart.z) * (float)step.z);

    // tDelta = abs(1/d)
    float3 tDelta = abs(1.0.xxx / max(abs(d), 1e-20.xxx)); // avoid div0

    // tMax as in listing: choose next boundary based on sign(d)
    float3 nextBoundary;
    nextBoundary.x = (d.x > 0.0) ? ceil(start.x) : floor(start.x);
    nextBoundary.y = (d.y > 0.0) ? ceil(start.y) : floor(start.y);
    nextBoundary.z = (d.z > 0.0) ? ceil(start.z) : floor(start.z);

    float3 tMax = (nextBoundary - start) / max(d, 1e-20.xxx);

    int3 voxel = (int3)floor(start);

    WriteOutput(voxel);
    GapDetection(voxel, flatAdvance, pAxis, normal, v0, v1, v2);

    // DDA loop
    [loop]
    while (stepsCount-- >= 0)
    {
        // find min(tMax)
        float minT = min(tMax.x, min(tMax.y, tMax.z));

        // TMaxStep = equal(tMax, minT)
        // Use epsilon compare for floats.
        const float eps = 1e-6;
        bool sx = abs(tMax.x - minT) <= eps;
        bool sy = abs(tMax.y - minT) <= eps;
        bool sz = abs(tMax.z - minT) <= eps;

        int3 TMaxStep = int3(sx ? 1 : 0, sy ? 1 : 0, sz ? 1 : 0);

        // Bias opposing the flat advance direction when multiple candidates exist:
        // Paper logic: if more than one axis ties, prefer voxel "behind" the advance direction.
        // In listing: mix(TMaxStep, notEqual(TMaxStep, flatAdvance), dot(TMaxStep,TMaxStep) > 1)
        int tieCount = TMaxStep.x + TMaxStep.y + TMaxStep.z;
        if (tieCount > 1)
        {
            // notEqual(TMaxStep, flatAdvance) per-component
            int3 ne;
            ne.x = (TMaxStep.x != flatAdvance.x) ? 1 : 0;
            ne.y = (TMaxStep.y != flatAdvance.y) ? 1 : 0;
            ne.z = (TMaxStep.z != flatAdvance.z) ? 1 : 0;
            TMaxStep = ne;
        }

        // Advance
        voxel += TMaxStep * step;

        // Update tMax only for stepped axes
        if (TMaxStep.x != 0) tMax.x += tDelta.x;
        if (TMaxStep.y != 0) tMax.y += tDelta.y;
        if (TMaxStep.z != 0) tMax.z += tDelta.z;

        WriteOutput(voxel);
        GapDetection(voxel, flatAdvance, pAxis, normal, v0, v1, v2);
    }
}

// Listing 1: VoxelizeTriangle overview (two halves + edge processing)
void VoxelizeTriangle(float3 v0, float3 v1, float3 v2)
{
    float3 normal = CalcNormal(v0, v1, v2);

    float3 pAxis = CalcPrimaryAxis(normal);

    SortVertices(v0, v1, v2, pAxis);

    int3 flatAdvance = CalcFlatAdvance(pAxis, normal);

    // Interior traversal from v2 up to v0 along primary axis (integer floors)
    int h = (int)floor(dot(v2, pAxis));
    int hMid = (int)floor(dot(v1, pAxis));
    int hTop = (int)floor(dot(v0, pAxis));

    // Bottom half: v2 -> v1 (scanlines between edges (v2,v0) and (v2,v1))
    [loop]
    while (h <= hMid)
    {
        float3 start = PointInEdge(v2, v0, pAxis, (float)h);
        float3 end   = PointInEdge(v2, v1, pAxis, (float)h);
        Scanline(start, end, flatAdvance, pAxis, normal, v0, v1, v2);
        h++;
    }

    // Top half: v1 -> v0 (scanlines between edges (v2,v0) and (v1,v0))
    [loop]
    while (h <= hTop)
    {
        float3 start = PointInEdge(v2, v0, pAxis, (float)h);
        float3 end   = PointInEdge(v1, v0, pAxis, (float)h);
        Scanline(start, end, flatAdvance, pAxis, normal, v0, v1, v2);
        h++;
    }

    // Edge voxelization (three edges)
    Scanline(v0, v1, flatAdvance, pAxis, normal, v0, v1, v2);
    Scanline(v1, v2, flatAdvance, pAxis, normal, v0, v1, v2);
    Scanline(v2, v0, flatAdvance, pAxis, normal, v0, v1, v2);
}

// -----------------------------
// Kernels
// -----------------------------

[numthreads(256, 1, 1)]
void ClearVolume(uint3 tid : SV_DispatchThreadID)
{
    uint idx = tid.x;
    uint total = (uint)(_Dim.x * _Dim.y * _Dim.z);
    if (idx >= total) return;
    _Voxels[idx] = 0;
}

[numthreads(128, 1, 1)]
void VoxelizeTriangles(uint3 tid : SV_DispatchThreadID)
{
    uint triId = tid.x;
    if (triId >= _TriCount) return;

    Tri t = _Tris[triId];

    // If you need transforms, apply _WorldToVoxel here before voxelization.
    VoxelizeTriangle(t.v0, t.v1, t.v2);
}

Notes (important when you wire this into Unity)

Dispatch sizes
- Clear: ceil(totalVoxels/256)
- Voxelize: ceil(triCount/128)
Atomic behavior
- Multiple triangles will race to set the same voxel; that’s fine because we only set to 1.
Performance
- The paper is “single compute pass” in GLSL. Unity will still be single pass, but the real perf depends on your triangle distribution and memory bandwidth, just like in their results. citeturn1view5turn1view6
Correctness details
- The paper’s gap detection is very specific: it checks a plane/edge intersection within the current floor and then does an inside-triangle test. citeturn1view3turn1view4
- If you feed non-watertight/very noisy triangles, you’ll still get a conservative-ish surface, but expect edge cases around degenerates.

If you want, I can also drop in a minimal C# setup snippet (buffers, dispatch, reading back/visualizing), and/or convert output to a packed bitset (8× memory reduction) while keeping atomics safe.

Below is a drop-in function you can add to the same Unity compute shader to remap world space → voxel space using a world-space AABB (min/max), exactly as you asked. I’ll also show where to call it in the pipeline so it cleanly integrates with the paper’s algorithm.

1) World → voxel remapping function (AABB-based)

Concept (matches voxelization literature)

We map a world position Pw into voxel grid coordinates such that:

bboxMin → (0,0,0)
bboxMax → (_Dim.x, _Dim.y, _Dim.z)
Coordinates are continuous (floats), not snapped
Voxel centers are implicitly at integer coordinates

Formula:

Pv = (Pw − bboxMin) / (bboxMax − bboxMin) * voxelResolution

2) HLSL function (Unity compute shader)

Add these uniforms:

float3 _BBoxMin;   // world-space bounding box minimum
float3 _BBoxMax;   // world-space bounding box maximum

Add this function:

// ----------------------------------
// World → Voxel space mapping
// ----------------------------------
float3 WorldToVoxel(float3 worldPos)
{
    float3 size = _BBoxMax - _BBoxMin;

    // Avoid division by zero on degenerate bounds
    size = max(size, 1e-8.xxx);

    // Normalize to [0,1]
    float3 normalized = (worldPos - _BBoxMin) / size;

    // Map to voxel grid coordinates
    // NOTE: This maps bboxMax exactly to Dim, not Dim-1.
    // floor() later will clamp naturally.
    return normalized * (float3)_Dim;
}

This matches what voxelizers typically do in CUDA/GLSL implementations and preserves sub-voxel precision for scanline traversal.

3) Where to use it (important)

🔴 DO NOT remap inside scanline logic

You must remap vertices once, before calling VoxelizeTriangle, otherwise you will distort the geometry.

Modify your kernel like this:

[numthreads(128, 1, 1)]
void VoxelizeTriangles(uint3 tid : SV_DispatchThreadID)
{
    uint triId = tid.x;
    if (triId >= _TriCount) return;

    Tri t = _Tris[triId];

    // World → voxel space
    float3 v0 = WorldToVoxel(t.v0);
    float3 v1 = WorldToVoxel(t.v1);
    float3 v2 = WorldToVoxel(t.v2);

    VoxelizeTriangle(v0, v1, v2);
}

That’s it — everything else in the paper stays unchanged, because:

Scanlines assume a uniform voxel grid
Gap detection relies on unit voxel edges
Amanatides & Woo traversal assumes axis-aligned voxels

4) Notes on correctness (very important)

A) Why we map to `[0, Dim]` instead of `[0, Dim−1]`

This avoids edge loss:

A triangle touching bboxMax still produces voxels
floor() later ensures valid voxel indices
This is the same convention used in conservative raster voxelizers

B) Clamping (optional)

If your meshes slightly exceed the bounding box:

float3 v = WorldToVoxel(worldPos);
v = clamp(v, 0.0.xxx, (float3)_Dim - 1e-4.xxx);

Avoid hard clamping unless you know you want cropping.

C) Precision

Use float, not half
Keep voxel space values continuous until floor()
This preserves the paper’s gap detection guarantees

5) Summary (what you now have)

You now have:

✅ Paper-accurate scanline + gap detection voxelizer
✅ Correct world → voxel remapping via AABB
✅ Compatible with Unity, DX12, Vulkan
✅ Safe for conservative surface voxelization

If you want next steps, I can:

Add solid voxel filling (GPU parity or flood fill)
Add bit-packed voxels (8×–32× memory reduction)
Add mesh batch voxelization (multiple AABBs / chunks)
Validate against paper failure cases (thin triangles, steep slopes)

Just tell me which direction you want to go.

Depends what you optimize for. For this voxelizer (lots of random writes, many collisions, “set this voxel to 1”), a RWStructuredBuffer is usually better. A RWTexture3D can be better mainly when you want easy sampling/filtering later or you’re doing mostly coherent writes.

Write performance and correctness

ComputeBuffer (RWStructuredBuffer)

Usually faster for random write patterns: it’s just a linear UAV store.
Atomics are straightforward (InterlockedExchange/Or/Add) and reliable on UAV buffers.
Easy to bit-pack (huge win for memory bandwidth), e.g. 1 bit/voxel with InterlockedOr on a uint word.

3D Texture (RWTexture3D / RenderTexture volume)

Random writes often go through texture tiling/swizzle layouts; that can be worse for scattered access.
Atomics support is more limited/quirky across APIs and Unity backends:
- HLSL has Interlocked* on RWTexture* in DX, but Unity/platform support can vary.
- Also depends on texture format (you’d want something like R32_UINT if you need atomics).
If you don’t need atomics (e.g., writing unique voxels), texture can be fine—but your voxelizer definitely has collisions.

Memory and bandwidth

Texture3D is typically heavier than a packed buffer if you store R8/R16 etc? It can be compact, but:
- Many platforms don’t like UAV writes to compressed/packed formats.
- Bitpacking in a buffer can beat almost any 3D texture format for a binary occupancy grid.
If you store R32_UINT 3D texture, it’s the same 4 bytes/voxel as a uint buffer, but often with more overhead in cache behavior.

Downstream usage (this is where Texture3D wins)

Pick Texture3D if you want:

Immediate use in shaders for raymarching, SDF-ish effects, 3D sampling, mipmapping (though mips for UAV-written volume usually need a separate pass), etc.
Cleaner integration with some rendering pipelines that expect a volume texture.

Pick ComputeBuffer if you want:

Further compute processing (flood fill, parity fill, morphology, meshing) where you’ll do more random reads/writes.
Bitpacked occupancy and fast atomics.
Fast CPU readback (often easier/cheaper than 3D texture readback paths).

Practical recommendation for your algorithm

Best overall pattern:

Voxelize into a buffer (optionally bitpacked with atomics).
If you need it as a volume for rendering, convert buffer → Texture3D in a second compute pass (coalesced writes, no collisions).

That gives you:

Fast/robust voxelization writes (buffer + atomics)
Easy sampling later (texture)

Unity-specific gotchas

Volume RenderTexture must be created with:
- dimension = TextureDimension.Tex3D
- enableRandomWrite = true
- a UAV-capable format (GraphicsFormat.R32_UInt if you want atomic-like integer ops; otherwise R8_UNorm/R16 etc for plain stores)
Even when it “works”, performance can vary a lot by GPU/driver.

Bottom line:
For this conservative scanline voxelizer: ComputeBuffer is typically better (faster + safer atomics + easier packing). Use a 3D texture only if the next stage strongly benefits from volume sampling, and ideally populate it in a separate conversion pass.