// 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
// -----------------------------
// -----------------------------
StructuredBuffer<Tri> _Tris;
RWStructuredBuffer<uint> _Voxels; // length = _Dim.x * _Dim.y * _Dim.z
// Optional: coordinate transform (if you feed object/world space)
// float4x4 _WorldToVoxel; // not used in core listing port
// -----------------------------
// -----------------------------
return (uint)(p.x + p.y * _Dim.x + p.z * _Dim.x * _Dim.y);
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));
// -----------------------------
// -----------------------------
// CalcPrimaryAxis: axis with smallest abs(normal) component (least dominant axis).
float3 CalcPrimaryAxis(float3 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);
// 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}
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);
// 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);
// Degenerate edge w.r.t pAxis: return vA (safe-ish)
if (abs(denom) < 1e-8) return vA;
float t = (h - a) / denom;
// 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 c0 = cross(e0, p - a);
float3 c1 = cross(e1, p - b);
float3 c2 = cross(e2, p - c);
// Allow small negative due to fp
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
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;
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)
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);
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)
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);
GapDetection(voxel, flatAdvance, pAxis, normal, v0, v1, v2);
while (stepsCount-- >= 0)
float minT = min(tMax.x, min(tMax.y, tMax.z));
// TMaxStep = equal(tMax, minT)
// Use epsilon compare for floats.
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;
// notEqual(TMaxStep, flatAdvance) per-component
ne.x = (TMaxStep.x != flatAdvance.x) ? 1 : 0;
ne.y = (TMaxStep.y != flatAdvance.y) ? 1 : 0;
ne.z = (TMaxStep.z != flatAdvance.z) ? 1 : 0;
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;
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))
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);
// Top half: v1 -> v0 (scanlines between edges (v2,v0) and (v1,v0))
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);
// 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);
// -----------------------------
// -----------------------------
void ClearVolume(uint3 tid : SV_DispatchThreadID)
uint total = (uint)(_Dim.x * _Dim.y * _Dim.z);
if (idx >= total) return;
void VoxelizeTriangles(uint3 tid : SV_DispatchThreadID)
if (triId >= _TriCount) return;
// If you need transforms, apply _WorldToVoxel here before voxelization.
VoxelizeTriangle(t.v0, t.v1, t.v2);