Differentiable Voronoi Ray Tracing

• Differentiable Voronoi ray tracing is a rendering method that partitions space into nearest-site cells and differentiates both within-cell quantities and moving boundaries for applications such as inverse rendering and novel-view synthesis.
• Its core pipeline traverses Delaunay-connected Voronoi cells along each ray, computes exact segment lengths and volume compositing, and propagates gradients through density, radiance, site positions, camera parameters, and ray geometry.
• Complete derivatives require boundary-integral estimators or softened assignments to capture ownership changes, visibility events, and topology transitions that ordinary autodiff misses, while recent systems report real-time performance exceeding 600 FPS in some settings.

Differentiable Voronoi ray tracing is a class of ray-based rendering methods in which space, a surface parameter domain, or a light-transport domain is partitioned into regions associated with nearest sites, and rendering is differentiated with respect to both the smooth quantities within each region and the moving boundaries between regions. Its defining difficulty is that nearest-site ownership, ray–cell traversal, visibility, and first-hit selection are discrete operations. Ordinary automatic differentiation differentiates only the currently active branch; a complete derivative must additionally account for boundary motion, changes in cell assignment, and topology events. Recent work spans three complementary approaches: explicit boundary-integral estimators for discontinuous rendering, piecewise-differentiable Voronoi geometry, and differentiable volumetric rendering through Voronoi-cell traversal, including real-time implementations such as Radiant Foam and VoroTracing (Li, 2019, Shumilin et al., 2023, Govindarajan et al., 3 Feb 2025, Taveira et al., 18 Aug 2026).

1. Mathematical formulation and derivative structure

Let sites depend on scene parameters $\Phi$ and define a Voronoi partition

$$V_i(\Phi)=\{x:d_i(x;\Phi)\leq d_j(x;\Phi)\ \forall j\}.$$

For Euclidean distance, $d_i(x;\Phi)=\|x-s_i(\Phi)\|$. A rendered quantity integrated over a domain $D$ can be expressed as

$$I(\Phi)=\sum_i\int_{D\cap V_i(\Phi)} f_i(x;\Phi)\,dx,$$

where $f_i$ is smooth while the active site remains fixed. Between sites $i$ and $j$, the cell boundary is defined implicitly by

$\alpha_{ij}(x;\Phi)=d_i(x;\Phi)-d_j(x;\Phi)=0.$

For squared Euclidean distance,

$\alpha_{ij}(x)=\|x-s_i\|^2-\|x-s_j\|^2,$

which expands to

$$V_i(\Phi)=\{x:d_i(x;\Phi)\leq d_j(x;\Phi)\ \forall j\}.$$0

Consequently, the bisector is planar and

$$V_i(\Phi)=\{x:d_i(x;\Phi)\leq d_j(x;\Phi)\ \forall j\}.$$1

The derivative of the integral has two components:

$$V_i(\Phi)=\{x:d_i(x;\Phi)\leq d_j(x;\Phi)\ \forall j\}.$$2

The first term is the smooth or interior derivative. It includes shading, radiance, materials, ray directions, textures, and other computations performed without changing nearest-site ownership. The second is the boundary term. It accounts for the motion of the Voronoi boundary and the jump in rendered contribution between neighboring cells.

This decomposition is the Voronoi analogue of the Heaviside/Dirac formulation used for visibility discontinuities in differentiable Monte Carlo ray tracing. A piecewise rendering function can be written locally as

$$V_i(\Phi)=\{x:d_i(x;\Phi)\leq d_j(x;\Phi)\ \forall j\}.$$3

where $$V_i(\Phi)=\{x:d_i(x;\Phi)\leq d_j(x;\Phi)\ \forall j\}.$$4 is the Heaviside function. Differentiating $$V_i(\Phi)=\{x:d_i(x;\Phi)\leq d_j(x;\Phi)\ \forall j\}.$$5 produces a Dirac delta,

$$V_i(\Phi)=\{x:d_i(x;\Phi)\leq d_j(x;\Phi)\ \forall j\}.$$6

which converts to a boundary integral. The same principle is equivalently expressed by Reynolds transport theorem: differentiating an integral over a moving cell yields an interior term plus a boundary flux term proportional to boundary velocity and the jump between neighboring-cell integrands (Li, 2019).

If $$V_i(\Phi)=\{x:d_i(x;\Phi)\leq d_j(x;\Phi)\ \forall j\}.$$7 on a common boundary, the boundary contribution vanishes even when the nearest-site assignment changes. If the sites produce different radiance, materials, densities, or geometric responses, the jump is nonzero and the boundary term must be represented.

2. Ray traversal through Voronoi cells

A ray is written as

$$V_i(\Phi)=\{x:d_i(x;\Phi)\leq d_j(x;\Phi)\ \forall j\}.$$8

where $$V_i(\Phi)=\{x:d_i(x;\Phi)\leq d_j(x;\Phi)\ \forall j\}.$$9 is the origin and $d_i(x;\Phi)=\|x-s_i(\Phi)\|$0 is a unit direction. Once the ray is inside the cell associated with site $d_i(x;\Phi)=\|x-s_i(\Phi)\|$1, only neighboring cells need to be examined. For a neighboring site $d_i(x;\Phi)=\|x-s_i(\Phi)\|$2, the shared bisector plane is

$d_i(x;\Phi)=\|x-s_i(\Phi)\|$3

Substitution of the ray equation gives the intersection parameter

$d_i(x;\Phi)=\|x-s_i(\Phi)\|$4

provided the denominator is nonzero. The next cell is selected by the smallest valid forward intersection among the current cell’s neighbors:

$d_i(x;\Phi)=\|x-s_i(\Phi)\|$5

If consecutive boundary intersections occur at $d_i(x;\Phi)=\|x-s_i(\Phi)\|$6 and $d_i(x;\Phi)=\|x-s_i(\Phi)\|$7, the segment length within the current cell is

$d_i(x;\Phi)=\|x-s_i(\Phi)\|$8

The traversal algorithm therefore consists of an initial nearest-site query followed by local face tests:

1. identify the initial cell;
2. enumerate its Delaunay-derived neighboring cells;
3. intersect the ray with candidate bisector planes;
4. reject invalid or non-forward intersections;
5. select the first valid exit;
6. integrate the current cell over the resulting interval;
7. move to the neighboring cell;
8. repeat until the ray exits the domain or transmittance becomes sufficiently small.

The computational advantage is locality: after initialization, traversal cost is governed primarily by the number of cells crossed and the degree of the visited cells rather than directly by the total number of sites. For a cell with $d_i(x;\Phi)=\|x-s_i(\Phi)\|$9 neighbors, the local exit operation costs $D$0. A ray crossing $D$1 cells has approximate traversal cost

$D$2

which is effectively $D$3 for bounded-degree meshes.

This traversal is implemented directly in “Radiant Foam: Real-Time Differentiable Ray Tracing,” which uses a 3D Voronoi tessellation as the rendered volumetric partition and a Delaunay structure as its dual adjacency representation (Govindarajan et al., 3 Feb 2025). “Differentiable Voronoi Ray Tracing Beyond Rasterization Speeds” develops the same principle into VoroTracing, emphasizing reduced traversal length, per-cell work, memory locality, and coherent GPU execution (Taveira et al., 18 Aug 2026).

3. Differentiable volume rendering

In Voronoi volumetric rendering, density and radiance are constant within each cell or are evaluated from cell-associated appearance parameters. For a segment of length $D$4 in cell $D$5, opacity is

$D$6

The transmittance before the segment is

$D$7

and its compositing weight is

$D$8

The rendered color is

$D$9

where $$I(\Phi)=\sum_i\int_{D\cap V_i(\Phi)} f_i(x;\Phi)\,dx,$$0 is a representative point and $$I(\Phi)=\sum_i\int_{D\cap V_i(\Phi)} f_i(x;\Phi)\,dx,$$1 is the cell radiance. The recurrence is

$$I(\Phi)=\sum_i\int_{D\cap V_i(\Phi)} f_i(x;\Phi)\,dx,$$2

For piecewise-constant density and radiance, the cellwise sum is the exact integral of the chosen representation rather than a numerical approximation obtained by sampling arbitrary points along the ray. Gradients propagate through density, segment lengths, alpha, transmittance, cell radiance, ray origins, ray directions, and site positions.

The opacity derivatives are

$$I(\Phi)=\sum_i\int_{D\cap V_i(\Phi)} f_i(x;\Phi)\,dx,$$3

and

$$I(\Phi)=\sum_i\int_{D\cap V_i(\Phi)} f_i(x;\Phi)\,dx,$$4

Because $$I(\Phi)=\sum_i\int_{D\cap V_i(\Phi)} f_i(x;\Phi)\,dx,$$5 depends on bisector intersections, site-position gradients reach the rendered image through the moving Voronoi geometry. A site displacement changes cell boundaries, ray exit times, segment lengths, surface-direction parameters, texture coordinates, and compositing weights.

Radiant Foam parameterizes density using a softplus activation with $$I(\Phi)=\sum_i\int_{D\cap V_i(\Phi)} f_i(x;\Phi)\,dx,$$6 and represents view-dependent color with degree-three spherical harmonics. VoroTracing instead optimizes an unconstrained density parameter $$I(\Phi)=\sum_i\int_{D\cap V_i(\Phi)} f_i(x;\Phi)\,dx,$$7 using

$$I(\Phi)=\sum_i\int_{D\cap V_i(\Phi)} f_i(x;\Phi)\,dx,$$8

Since

$$I(\Phi)=\sum_i\int_{D\cap V_i(\Phi)} f_i(x;\Phi)\,dx,$$9

the density gradient becomes

$f_i$0

The segment length cancels, making the density optimization scale-invariant with respect to cell size. In the reported ablation, exponential density improved average PSNR from $f_i$1 to $f_i$2 dB, reduced cells per ray from $f_i$3 to $f_i$4, and increased speed from $f_i$5 to $f_i$6 FPS (Taveira et al., 18 Aug 2026).

VoroTracing uses two $f_i$7 RGB octahedral textures per cell: a view-independent surface texture and a view-dependent residual. Given a direction $f_i$8 from the site to a cell-boundary point,

$f_i$9

and the direction is unfolded onto a square through octahedral mapping. The final color is obtained by applying a componentwise sigmoid to the sum of the surface-dependent and view-dependent texture lookups. Bilinear interpolation distributes gradients to the four neighboring texels of each map.

The representation contrasts with spherical harmonics. Degree-three spherical harmonics require $i$0 RGB coefficient values per cell, whereas two bilinear octahedral lookups load $i$1 values in the nominal case. The textures also encode spatial variation across a cell’s surface footprint, allowing larger cells to carry more appearance detail without necessarily increasing traversal depth.

4. Voronoi construction, autodifferentiation, and topology

A differentiable 2D Voronoi construction can be obtained by treating the Delaunay triangulation as discrete external information and differentiating the resulting geometric operations. For a Delaunay triangle with sites $i$2, $i$3, and $i$4, the corresponding Voronoi vertex is the circumcenter. Its coordinates are rational functions of the site coordinates, with denominator

$i$5

The circumcenter is differentiable wherever $i$6. Each interior Delaunay edge corresponds to a Voronoi edge joining the circumcenters of its two incident triangles. A Voronoi cell is formed by ordering the circumcenters and, for unbounded cells, adding differentiable ghost points. Cell areas can be computed with the Shoelace formula, and bounded diagrams can be produced by differentiable polygon clipping (Shumilin et al., 2023).

The resulting computational structure is piecewise differentiable:

$i$7

The Delaunay triangulation itself is not included in the autodiff graph. It is recomputed periodically or updated incrementally, while gradients are propagated through the current geometric graph. This supplies local derivatives for fixed topology but not derivatives through Delaunay flips, cocircular configurations, or changes in which cells are unbounded.

Radiant Foam likewise maintains Delaunay-derived adjacency while rendering Voronoi cells. Its argument is that a Voronoi face disappears with zero area during an adjacency transition, so the piecewise-constant volumetric field and ray-segment lengths remain continuous under ordinary site motion except at degenerate configurations. The exact derivative at the topology event is nevertheless not generally defined (Govindarajan et al., 3 Feb 2025).

Important degeneracies include nearly collinear or coplanar site configurations, coincident sites, zero-area Voronoi faces, multiple-site ties, ray–face parallelism, changes in the first intersected face, and simultaneous boundary events. Circumcenter derivatives become ill-conditioned when $i$8 approaches zero. The hard minimum selecting the next face and the hard argmin selecting the nearest site are nondifferentiable at ties. A piecewise-autodiff implementation can accept these events as measure-zero configurations, but numerical robustness may require explicit degeneracy detection, stable predicates, regularization, symbolic perturbation, or event-aware graph rebuilding.

5. Boundary-aware and softened differentiation

The piecewise-differentiable Voronoi traversal used by Radiant Foam and VoroTracing differentiates ray–bisector intersections conditional on a fixed cell adjacency and active face. It does not, by itself, provide a globally smooth derivative through every change in cell identity. Two broad strategies address this limitation.

Boundary-integral estimators

The boundary-integral approach explicitly samples moving discontinuities. For an edge or boundary $i$9, a Monte Carlo estimator has the form

$j$0

The two-sided contributions are evaluated from slightly offset samples on either side of the boundary, typically with an offset of approximately $j$1 and a shared random-number sequence. This correlated sampling reduces variance. The method is not finite differencing: the scene parameters are not perturbed and rerendered; instead, the jump in contribution across a fixed geometric boundary is estimated directly.

For Voronoi ray tracing, triangle edges are replaced by active bisectors or higher-dimensional assignment boundaries. Candidate neighboring site pairs must be enumerated efficiently, and the boundary measure and Jacobian must be adapted to the selected metric. In screen space, a boundary is generally a curve; in ray space, it may be a hypersurface. The boundary contribution is proportional to the jump in complete rendered or transport contribution, not merely to the difference in site labels.

Soft partitions and continuation

A hard nearest-site assignment can be replaced by soft weights,

$j$2

or by a soft minimum,

$j$3

Soft assignments provide gradients across cell boundaries but modify the exact Voronoi model near those boundaries. The same trade-off appears in differentiable optical and radio ray tracing. Algorithmic differentiable non-sequential ray tracing uses automatic differentiation through ray states, refraction, detector intersections, and smooth reconstruction filters, while leaving discrete intersection and branch decisions only partially treated (Koning et al., 2023). Differentiable radio ray tracing differentiates physical field computation along fixed paths but does not differentiate path generation, visibility, geometry, or path topology (Hoydis et al., 2023).

RayLoc addresses sparse gradients and local minima by convolving its localization loss with a Gaussian over candidate positions and using coarse-to-fine variance decay. Without Gaussian smoothing, the device-free median error exceeds $j$4 m, compared with approximately $j$5 m for RayLoc. This strategy is relevant as a continuation method for Voronoi ownership and boundary events, but it is an optimization approximation rather than an exact boundary derivative (Han et al., 21 Jan 2025).

6. Systems, applications, and performance

Voronoi ray tracing has applications in differentiable novel-view synthesis, inverse rendering, optical design, radio-environment calibration, radiative transfer, and scene localization. The physical quantities that can be differentiated depend on the renderer. They may include site positions, camera parameters, ray origins and directions, densities, radiance textures, material properties, antenna patterns, propagation delays, reflection and diffraction coefficients, and task-specific losses.

For physically based rendering, the full differentiable pipeline separates smooth path derivatives from visibility or ownership derivatives. The general-purpose differentiable ray tracer in “Differentiable Visual Computing” supports camera pose, mesh geometry, lights, materials, textures, and arbitrary light-transport parameters, with a reported overhead of approximately $j$6–$j$7 over rendering an image alone. Its experiments report agreement with central finite differences within $j$8 relative to the $j$9 norm at $\alpha_{ij}(x;\Phi)=d_i(x;\Phi)-d_j(x;\Phi)=0.$0 resolution and demonstrate optimization involving primary visibility, shadows, glossy reflection, global illumination, and camera parameters (Li, 2019).

For volumetric Voronoi rendering, Radiant Foam reports $\alpha_{ij}(x;\Phi)=d_i(x;\Phi)-d_j(x;\Phi)=0.$1 FPS on Mip-NeRF 360 and $\alpha_{ij}(x;\Phi)=d_i(x;\Phi)-d_j(x;\Phi)=0.$2 FPS on Deep Blending, with average values of $\alpha_{ij}(x;\Phi)=d_i(x;\Phi)-d_j(x;\Phi)=0.$3 PSNR, $\alpha_{ij}(x;\Phi)=d_i(x;\Phi)-d_j(x;\Phi)=0.$4 SSIM, and $\alpha_{ij}(x;\Phi)=d_i(x;\Phi)-d_j(x;\Phi)=0.$5 LPIPS on Mip-NeRF 360. Its core engineering choices are local convex-cell traversal, exact piecewise-constant volume compositing, GPU execution without RT cores or OptiX, and Delaunay-derived adjacency (Govindarajan et al., 3 Feb 2025).

VoroTracing co-designs representation and execution. It concentrates opacity near surfaces, uses exponential density, employs compact octahedral textures, avoids pruning and densification with a fixed budget of $\alpha_{ij}(x;\Phi)=d_i(x;\Phi)-d_j(x;\Phi)=0.$6 million sites, and applies Morton ordering, $\alpha_{ij}(x;\Phi)=d_i(x;\Phi)-d_j(x;\Phi)=0.$7 warp-coherent ray tiles, aligned four-channel texture loads, half-precision attributes, and low-contribution cell skipping. On Mip-NeRF 360, it reports $\alpha_{ij}(x;\Phi)=d_i(x;\Phi)-d_j(x;\Phi)=0.$8 FPS on an RTX 5090, compared with $\alpha_{ij}(x;\Phi)=d_i(x;\Phi)-d_j(x;\Phi)=0.$9 FPS for Radiant Foam, $\alpha_{ij}(x)=\|x-s_i\|^2-\|x-s_j\|^2,$0 FPS for ray-traced PowerFoam, and $\alpha_{ij}(x)=\|x-s_i\|^2-\|x-s_j\|^2,$1 FPS for 3D Gaussian Splatting. Its average quality is $\alpha_{ij}(x)=\|x-s_i\|^2-\|x-s_j\|^2,$2 PSNR, $\alpha_{ij}(x)=\|x-s_i\|^2-\|x-s_j\|^2,$3 SSIM, and $\alpha_{ij}(x)=\|x-s_i\|^2-\|x-s_j\|^2,$4 LPIPS. The reported average traversal length is $\alpha_{ij}(x)=\|x-s_i\|^2-\|x-s_j\|^2,$5 cells per ray, compared with $\alpha_{ij}(x)=\|x-s_i\|^2-\|x-s_j\|^2,$6 for Radiant Foam (Taveira et al., 18 Aug 2026).

Its ray formulation naturally supports non-pinhole effects by changing ray generation rather than redesigning a rasterizer. Fisheye projection, rolling shutter, motion blur, and depth of field can be expressed through altered ray origins, directions, camera poses, or multisampling. The reported $\alpha_{ij}(x)=\|x-s_i\|^2-\|x-s_j\|^2,$7 FPS applies to standard pinhole rendering; multisampled effects reduce frame rate in proportion to the number of rays per pixel.

A complete differentiable Voronoi ray tracer must nevertheless address several unresolved or conditionally solved issues:

• Nearest-site ownership: hard argmin operations are differentiable only away from ties.
• First-face selection: the minimum over candidate intersections changes discontinuously when two faces exchange order.
• Topology changes: Delaunay flips and Voronoi adjacency changes require graph updates and do not receive ordinary reverse-mode derivatives.
• Degeneracy: coincident sites, zero-area faces, nearly singular circumcenters, and parallel ray–face configurations can produce unstable Jacobians.
• Visibility and path birth/death: exact geometric events require boundary-integral estimators, soft visibility, or other measure-aware treatments.
• Dynamic structure: pruning, densification, site insertion, and deletion introduce discrete or trans-dimensional changes beyond ordinary continuous differentiation.
• Physical interfaces: reflection, refraction, scattering, and global illumination require additional interface normals, material models, and secondary-ray estimators.

Ray-trax illustrates a related GPU architecture based on vectorized rays, fixed-size loops, JAX compilation, multi-GPU sharding, and differentiable emission–absorption updates, but it operates on a regular Cartesian grid with fixed spatial steps rather than Voronoi cells (Branca et al., 12 Nov 2025). Power Foam extends foam-based ray tracing toward bounded power diagrams and oriented surface formulations, but the supplied source does not provide technical details sufficient to establish its algorithm or results (Govindarajan et al., 27 Apr 2026).

The central methodological distinction is therefore between differentiating a smooth rendering computation conditional on a fixed spatial partition and differentiating the partition itself. Voronoi ray tracing provides explicit cell geometry, local traversal, and differentiable segment lengths; boundary-aware methods provide the missing contribution when ownership changes alter the integrated rendering function. A system that combines both—piecewise autodiff for fixed topology, explicit bisector derivatives for moving boundaries, and robust treatment of topology and visibility events—constitutes the technically complete form of differentiable Voronoi ray tracing.