voxelize

GPU voxelization for triangle meshes and point clouds, loadable through kernels. Meshes rasterize into bit-packed occupancy grids along two routes: a conservative surface (every cell a triangle touches) and a parity-filled solid (every cell whose center the mesh contains). Point clouds bin into occupancy grids or pooled per-voxel features. All grid outputs are bitwise deterministic: the same inputs produce the same bytes on every card, in any triangle or point order.

Usage

from kernels import get_kernel

vox = get_kernel("phanerozoic/voxelize", version=1, trust_remote_code=True)

grid = vox.voxelize_mesh(verts, faces, resolution=512, mode="solid")
grid.occupancy                            # bool [DZ, DY, DX]
grid.origin, grid.voxel_size, grid.dims   # world transform

occ = vox.voxelize_points(points, voxel_size=0.1)
pv = vox.pool_points(points, feats, voxel_size=0.1, reduce="mean")
pv.coords, pv.features, pv.counts         # int32 [M,3], fp32 [M,C], int32 [M]

version selects the release branch; trust_remote_code is required by kernels for publishers without the trusted-publisher mark.

vertices is float32 [V, 3] and faces integer [F, 3]; both move to the GPU if needed. mode is 'surface', 'solid', or 'conservative' (their union). With only resolution the grid is fitted around the mesh with pad empty voxels on each side; pass voxel_size (and optionally origin and dims) for caller-controlled grids. packed=True returns the raw int32 bit grid [ceil(DZ/32), DY, DX] with voxel (x, y, z) at word [z//32, y, x], bit z%32.

API

How it works

The surface pass tests each triangle against every cell of its grid-space bounding box with the Akenine-Möller separating-axis test and sets occupied bits with atomicOr. Triangles are bucketed by bounding-box cell count and processed by a single thread, a warp, or a whole block, so a handful of large polygons cannot serialize a launch full of small ones.

The solid pass computes crossing parity per grid column. Each triangle XORs one bit at the first voxel center above its intersection with every column inside its 2D projection, and a column scan then converts crossing flags to occupancy with a prefix-XOR (five shift-XOR steps per 32 voxels). The projected inside test runs in snapped int64 fixed point (1/256 voxel) with a top-left fill rule; integer edge functions negate exactly between the two triangles sharing an edge, so every boundary column is claimed by exactly one of them and watertight parity cannot leak.

The point pass is one atomicOr per point for occupancy. Feature pooling (dynamic voxelization, no per-voxel point cap) stable-sorts points by voxel key and reduces each voxel's features sequentially in that fixed order, one thread per voxel-channel pair, with mean, max, or sum.

Occupancy is written only through integer OR/XOR atomics, which commute exactly, and every floating-point step in the geometric predicates is an explicit round-per-op intrinsic, so ptxas cannot contract expressions differently per target. Grids are therefore invariant to triangle order, point order, launch configuration, and GPU architecture, and pooled features are bit-identical for a fixed input order.

Correctness

An axis-aligned box placed off voxel boundaries voxelizes to its exact analytic counts (58³ solid and 58³ − 56³ surface cells at the test placement). An icosphere's parity fill matches its divergence-theorem volume to 2 parts in 10⁴ at 128³ and agrees with trimesh's ray-cast contains on every sampled voxel; a torus lands within a part in 10³. Pooled features match a sequential fp32 reference bitwise in mean, max, and sum modes. Shuffling faces or points changes no output bit.

The test suite pins SHA-256 digests of five canonical grids and two pooled outputs, built from integer-hash points and an IEEE-basic-ops icosphere so the inputs are identical on every host. Digest equality on another machine certifies bit-identical voxelization across CUDA architectures, torch releases, and CUDA toolchains. The digests, minted on an RTX 6000 Ada (sm89), reproduce byte-for-byte on L4 (Ada sm89), H200 (Hopper sm90), and RTX PRO 6000 (Blackwell sm120), and on a GeForce RTX 3070 Ti laptop (Ampere sm86): four architecture generations, from a consumer laptop GPU to a datacenter card.

Measured

RTX 6000 Ada (sm89), event-timed medians of the voxelization kernels (mesh already on device, grid pre-allocated, file output excluded, to match how the compared tools report their own kernel time).

Surface voxelization against cuda_voxelizer v0.6, the standard fast GPU voxelizer, which is also a deterministic integer-atomic bit-grid method, so this is a pure speed comparison:

The gap widens with triangle count: at 1.31M triangles the kernel sustains 7.7 billion triangles per second (11.5 for solid). Below ~50K triangles both tools are launch-overhead-bound in the tens of microseconds and the ratio is not meaningful. Solid margins are larger still, but part of that is cuda_voxelizer's solid fill being slower than its surface path. trimesh 4.11 (CPU) voxelizes the same meshes three to four orders of magnitude slower.

The same kernel on a GeForce RTX 3070 Ti laptop (Ampere, 8 GB) voxelizes the 1.31M-triangle mesh into a 512³ grid in 0.74 ms and bins 5M points at 8 billion per second: about 4x slower than the Ada, tracking its core count, yet still faster than cuda_voxelizer on the workstation Ada and bit-identical to it.

Point pooling into voxels (mean, 32 feature channels), against the unique + index_add_ recipe that voxel feature encoders use:

index_add_ uses CUDA float atomics, so the torch default returns different bytes on every run: three identical calls produced three different outputs. Making it reproducible with torch.use_deterministic_algorithms(True) runs about 2x slower than voxelize. voxelize is faster than the nondeterministic default and matches a sequential fp32 reference to the bit. Occupancy binning runs at 12 to 13 billion points per second.

Requirements and limits

Float32 geometry, grids up to 8192 voxels per axis, one geometry per call, single device. Solid and conservative modes assume a watertight mesh (crossing parity of an open mesh is undefined); surface mode takes any triangle soup, and degenerate triangles voxelize conservatively as their bounding cells. Mean and sum pooling are deterministic for a fixed input point order; max pooling is order-invariant outright.

License

Apache-2.0.