Differentiable 2D Gaussian Ray Tracing

• The paper introduces a differentiable method that leverages 2D Gaussian disk primitives to enable efficient, ray-centric rendering, enhancing quality over traditional 3D Gaussian splatting.
• It employs BVH-accelerated ray traversal with compact proxy geometries and chunked intersection processing to accurately model sensor attributes and indirect lighting.
• The technique integrates front-to-back alpha-compositing and differentiable accumulation, enabling end-to-end optimization of geometric and radiometric parameters for LiDAR and inverse rendering tasks.

Differentiable 2D Gaussian ray tracing denotes a class of Gaussian rendering methods in which rays are traced against disk-like 2D Gaussian primitives in 3D, the resulting interactions are accumulated in front-to-back order, and the entire process remains differentiable with respect to geometric and radiometric or sensor-response parameters. In contrast to conventional 3D Gaussian splatting, which is centered on projecting Gaussians into screen space and rasterizing them, this formulation is explicitly ray-centric: it uses BVH-accelerated traversal, ray–primitive intersection, and transmittance or alpha-compositing along the ray. In late-2024 work, the formulation appears in two distinct but closely related settings: dynamic LiDAR re-simulation in LiDAR-RT (Zhou et al., 2024) and inverse rendering with inter-reflection in IRGS (Gu et al., 2024).

1. Conceptual lineage and relation to Gaussian rendering

Gaussian reconstruction kernels were proposed by Westover (1990) and were studied by the computer graphics community in the 1990s as an alternative geometric representation to meshes and point clouds. A modern differentiable rendering precursor is VoGE, which uses explicit 3D Gaussian ellipsoids and renders them by tracing rays through volumetric density, rather than by rasterizing triangles or points (Wang et al., 2022). In VoGE, each ellipsoidal Gaussian induces a 1D Gaussian along a ray, and the renderer approximates the contribution of primitive $k$ as

$$W_{p,k}=T(l_k)e^{q_k}, \qquad C(\mathbf{r})\approx \sum_{k=1}^K W_{p,k}\mathbf{c}_k,$$

with transmittance expressed in closed form through the error function. This established a ray-traced, transmittance-based alternative to rasterization for Gaussian primitives.

Differentiable 2D Gaussian ray tracing departs from that volumetric 3D-Gaussian view in a specific way. A common misconception is to treat “2D Gaussian” as referring only to projected screen-space Gaussians used during rasterization. In the relevant late-2024 formulations, the term also refers to a surface-like primitive in 3D: LiDAR-RT treats 2D Gaussians as planar disks for ray tracing, and IRGS defines each 2D Gaussian as a disk-like element with a supporting plane and an analytic ray–splat intersection (Zhou et al., 2024). IRGS makes the motivation explicit: 3DGS is rasterization-centric, the rendering equation requires visibility and indirect radiance along arbitrary incident directions, and direct ray tracing of pretrained 3D Gaussian primitives is not ideal because the ray–splat intersection is ambiguous and can cause noticeable rendering degradation (Gu et al., 2024).

This shift is significant because it changes the role of Gaussians from purely image-plane splats to explicit ray-interaction primitives. LiDAR-RT further argues that rasterization is ill-suited for the cylindrical LiDAR range image and for LiDAR’s active sensing process, so ray tracing is adopted as the physically accurate core (Zhou et al., 2024).

2. Primitive parameterization and geometric proxies

IRGS parameterizes each 2D Gaussian by a center $\boldsymbol{\mu}\in\mathbb{R}^3$, opacity $o\in[0,1]$, two tangent vectors $\boldsymbol{t_u},\boldsymbol{t_v}\in\mathbb{R}^3$, and a scaling vector $\boldsymbol{s}=(s_u,s_v)\in\mathbb{R}^2$ (Gu et al., 2024). These parameters define an oriented disk-like primitive in 3D. The associated normal is $$\boldsymbol{n}=\boldsymbol{t_u}\times \boldsymbol{t_v}$$, which gives the primitive a well-defined supporting plane.

LiDAR-RT begins from the standard 3D Gaussian splatting primitive

$$\mathcal{G}_{i}(\mathbf{x}) = \exp\left(-\frac{1}{2}(\mathbf{x}-\mu_i)^\top \Sigma_i^{-1}(\mathbf{x}-\mu_i)\right), \qquad \Sigma_i = R_i S_i S_i^\top R_i^\top,$$

with mean position $\mu_i$, covariance $\Sigma_i$, opacity $$W_{p,k}=T(l_k)e^{q_k}, \qquad C(\mathbf{r})\approx \sum_{k=1}^K W_{p,k}\mathbf{c}_k,$$0, and spherical harmonics coefficients $$W_{p,k}=T(l_k)e^{q_k}, \qquad C(\mathbf{r})\approx \sum_{k=1}^K W_{p,k}\mathbf{c}_k,$$1. It then extends each Gaussian with LiDAR-specific learnable attributes $$W_{p,k}=T(l_k)e^{q_k}, \qquad C(\mathbf{r})\approx \sum_{k=1}^K W_{p,k}\mathbf{c}_k,$$2, where $$W_{p,k}=T(l_k)e^{q_k}, \qquad C(\mathbf{r})\approx \sum_{k=1}^K W_{p,k}\mathbf{c}_k,$$3 is reflection intensity and $$W_{p,k}=T(l_k)e^{q_k}, \qquad C(\mathbf{r})\approx \sum_{k=1}^K W_{p,k}\mathbf{c}_k,$$4 is ray-drop probability; both are modeled with spherical harmonics because intensity and drop probability depend on viewing direction (Zhou et al., 2024).

Efficient ray tracing requires proxy geometry. LiDAR-RT does not intersect rays against full analytic Gaussian volumes directly. Instead, it prefers 2D Gaussians as the ray-traced primitive, treats them as planar disks, and uses a pair of co-planar triangles as the simplest proxy geometry for BVH construction and intersection testing. The paper states that this is the most efficient option, because it reduces mesh complexity, tightly wraps the Gaussian, avoids the “maximum response” approximation used in some particle ray tracing methods, and makes the ray intersection sample locations coincide with actual ray hits (Zhou et al., 2024).

IRGS adopts a different proxy strategy. Each 2D Gaussian is enclosed by an adaptive icosahedron mesh sized to cover the primitive’s influence down to a minimum opacity threshold $$W_{p,k}=T(l_k)e^{q_k}, \qquad C(\mathbf{r})\approx \sum_{k=1}^K W_{p,k}\mathbf{c}_k,$$5: $$W_{p,k}=T(l_k)e^{q_k}, \qquad C(\mathbf{r})\approx \sum_{k=1}^K W_{p,k}\mathbf{c}_k,$$6 Each icosahedron has 20 triangular faces, so a scene with $$W_{p,k}=T(l_k)e^{q_k}, \qquad C(\mathbf{r})\approx \sum_{k=1}^K W_{p,k}\mathbf{c}_k,$$7 Gaussians yields $$W_{p,k}=T(l_k)e^{q_k}, \qquad C(\mathbf{r})\approx \sum_{k=1}^K W_{p,k}\mathbf{c}_k,$$8 triangles before BVH construction (Gu et al., 2024).

The coexistence of triangle-pair and icosahedral proxies indicates that differentiable 2D Gaussian ray tracing is not tied to a single geometric surrogate. What is shared is the use of compact bounding geometry to enable hardware-accelerated ray traversal over surface-like Gaussian elements.

3. Ray traversal, ordering, and intersection evaluation

In LiDAR-RT, rays are batched as

$$W_{p,k}=T(l_k)e^{q_k}, \qquad C(\mathbf{r})\approx \sum_{k=1}^K W_{p,k}\mathbf{c}_k,$$9

where $\boldsymbol{\mu}\in\mathbb{R}^3$0 are ray origins and $\boldsymbol{\mu}\in\mathbb{R}^3$1 are normalized ray directions. Ray generation is launched through NVIDIA OptiX using optixLaunch. The ray-generation program casts the rays against the BVH, the any-hit program records intersections, and a sorted buffer maintains the front-to-back order of intersections. To avoid the cost of storing and sorting all hits globally, intersections are processed in chunks; the paper uses chunk size 16. When a chunk fills, the renderer extracts indices $\boldsymbol{\mu}\in\mathbb{R}^3$2 of intersected primitives and depth values $\boldsymbol{\mu}\in\mathbb{R}^3$3, evaluates Gaussian responses at the sample points, accumulates the LiDAR attributes front-to-back, and advances until all intersections have been processed or the transmittance drops below $\boldsymbol{\mu}\in\mathbb{R}^3$4 (Zhou et al., 2024).

IRGS uses OptiX as well, but its emphasis is the exact ray–splat intersection for 2D Gaussians. Given ray origin $\boldsymbol{\mu}\in\mathbb{R}^3$5, direction $\boldsymbol{\mu}\in\mathbb{R}^3$6, Gaussian center $\boldsymbol{\mu}\in\mathbb{R}^3$7, and normal $\boldsymbol{\mu}\in\mathbb{R}^3$8, the intersection point is

$\boldsymbol{\mu}\in\mathbb{R}^3$9

The local coordinates in the Gaussian basis are

$o\in[0,1]$0

and the Gaussian response at the hit point is

$o\in[0,1]$1

A $o\in[0,1]$2-buffer per-ray sorting algorithm is then used to obtain the exact front-to-back ordering of intersected Gaussians along each ray (Gu et al., 2024).

These two formulations share a common operational structure: BVH traversal identifies candidate Gaussian proxies, per-ray ordering is established explicitly, and Gaussian responses are evaluated at ray-consistent sample points. The principal difference is task-driven. LiDAR-RT traces sensor rays to render range, intensity, and ray-drop; IRGS traces arbitrary incident-light rays to recover visibility and indirect radiance.

4. Accumulation, visibility, and differentiability

The accumulation model derives from the alpha-compositing structure used in Gaussian splatting preliminaries. For projected 2D Gaussians in screen space,

$o\in[0,1]$3

and the pixel value is

$o\in[0,1]$4

LiDAR-RT reuses the same compositing idea, but the rendered attribute is not RGB. Instead, $o\in[0,1]$5 is replaced by LiDAR-specific quantities such as intensity $o\in[0,1]$6 and ray-drop probability $o\in[0,1]$7, while depth or range is accumulated along the same front-to-back traversal (Zhou et al., 2024).

Its backward pass is organized to match the forward ordering. Rather than storing the full global blending order, LiDAR-RT differentiates in the same front-to-back order as the forward pass, using

$o\in[0,1]$8

where $o\in[0,1]$9 is the final rendered value and $\boldsymbol{t_u},\boldsymbol{t_v}\in\mathbb{R}^3$0 is the accumulated value up to the $\boldsymbol{t_u},\boldsymbol{t_v}\in\mathbb{R}^3$1-th Gaussian. The implementation launches a new optixLaunch, re-casts the same rays, re-obtains the sorted intersections, computes gradients with the front-to-back formula, and accumulates gradients in global buffers via atomic operations (Zhou et al., 2024). This is presented as crucial for differentiable rendering efficiency because it avoids the memory-heavy backward ordering used in camera-space 3DGS.

IRGS uses alpha blending in an analogous front-to-back form: $\boldsymbol{t_u},\boldsymbol{t_v}\in\mathbb{R}^3$2 In tracing mode,

$\boldsymbol{t_u},\boldsymbol{t_v}\in\mathbb{R}^3$3

and, for a sampled incident direction $\boldsymbol{t_u},\boldsymbol{t_v}\in\mathbb{R}^3$4 from surface point $\boldsymbol{t_u},\boldsymbol{t_v}\in\mathbb{R}^3$5,

$\boldsymbol{t_u},\boldsymbol{t_v}\in\mathbb{R}^3$6

The accumulated opacity therefore gives $\boldsymbol{t_u},\boldsymbol{t_v}\in\mathbb{R}^3$7, while the accumulated color gives indirect radiance $\boldsymbol{t_u},\boldsymbol{t_v}\in\mathbb{R}^3$8 (Gu et al., 2024). Because the outputs depend on Gaussian parameters, per-Gaussian view-dependent color, sorted ordering, and alpha blending, the method is fully differentiable and enables optimization of indirect light through backpropagation.

A central interpretive point follows directly from these formulations. Differentiable 2D Gaussian ray tracing is not merely a visibility oracle; it is a differentiable mechanism for aggregating task-specific attributes along a ray. In LiDAR-RT those attributes are depth, intensity, and ray-drop. In IRGS they are visibility and indirect radiance.

5. Task-specific instantiations

In LiDAR-RT, differentiable 2D Gaussian ray tracing is embedded in a dynamic-scene LiDAR re-simulation pipeline. The scene is decomposed into a static background and multiple foreground vehicles, and each object has its own Gaussian set. For dynamic objects, Gaussian parameters are defined in object-local coordinates and mapped into world coordinates by a tracked pose sequence $\boldsymbol{t_u},\boldsymbol{t_v}\in\mathbb{R}^3$9: $\boldsymbol{s}=(s_u,s_v)\in\mathbb{R}^2$0 This functions as the scene-graph-like update mechanism. Implementation details include background point-cloud fusion across frames, KNN-based normal estimation for Gaussian orientation initialization, voxel downsampling with voxel size $\boldsymbol{s}=(s_u,s_v)\in\mathbb{R}^2$1, and augmentation of object models with fewer than $\boldsymbol{s}=(s_u,s_v)\in\mathbb{R}^2$2 points by random sampling inside the 3D bounding box until they reach $\boldsymbol{s}=(s_u,s_v)\in\mathbb{R}^2$3 points (Zhou et al., 2024).

LiDAR output is rendered as a range image. For a point $\boldsymbol{s}=(s_u,s_v)\in\mathbb{R}^2$4,

$\boldsymbol{s}=(s_u,s_v)\in\mathbb{R}^2$5

The paper emphasizes, however, that the range image is the final representation rather than the core physical model: the physically accurate component is the ray-traced active-sensing simulation (Zhou et al., 2024). Training uses the total loss

$\boldsymbol{s}=(s_u,s_v)\in\mathbb{R}^2$6

with $\boldsymbol{s}=(s_u,s_v)\in\mathbb{R}^2$7, $\boldsymbol{s}=(s_u,s_v)\in\mathbb{R}^2$8, $\boldsymbol{s}=(s_u,s_v)\in\mathbb{R}^2$9, and $$\boldsymbol{n}=\boldsymbol{t_u}\times \boldsymbol{t_v}$$0. The loss terms are depth/range $$\boldsymbol{n}=\boldsymbol{t_u}\times \boldsymbol{t_v}$$1, intensity $$\boldsymbol{n}=\boldsymbol{t_u}\times \boldsymbol{t_v}$$2, ray-drop binary cross-entropy, and Chamfer Distance between rendered and ground-truth point clouds. After Gaussian optimization, a U-Net refines sensor-level ray drop using ray origin, ray direction, and predicted depth, intensity, and ray-drop (Zhou et al., 2024).

IRGS uses differentiable 2D Gaussian ray tracing inside the full rendering equation for inverse rendering: $$\boldsymbol{n}=\boldsymbol{t_u}\times \boldsymbol{t_v}$$3 The BRDF is a simplified Disney-style dielectric BRDF with diffuse albedo $$\boldsymbol{n}=\boldsymbol{t_u}\times \boldsymbol{t_v}$$4 and roughness $$\boldsymbol{n}=\boldsymbol{t_u}\times \boldsymbol{t_v}$$5, decomposed into $$\boldsymbol{n}=\boldsymbol{t_u}\times \boldsymbol{t_v}$$6, with $$\boldsymbol{n}=\boldsymbol{t_u}\times \boldsymbol{t_v}$$7. Incident radiance is decomposed as

$$\boldsymbol{n}=\boldsymbol{t_u}\times \boldsymbol{t_v}$$8

where direct lighting comes from an environment cubemap, and both visibility $$\boldsymbol{n}=\boldsymbol{t_u}\times \boldsymbol{t_v}$$9 and indirect radiance $$\mathcal{G}_{i}(\mathbf{x}) = \exp\left(-\frac{1}{2}(\mathbf{x}-\mu_i)^\top \Sigma_i^{-1}(\mathbf{x}-\mu_i)\right), \qquad \Sigma_i = R_i S_i S_i^\top R_i^\top,$$0 are queried on the fly by 2D Gaussian ray tracing (Gu et al., 2024).

IRGS uses a two-stage pipeline. Stage I pretrains a standard 2D Gaussian splatting model and optimizes geometry using RGB reconstruction, normal consistency, depth distortion, edge-aware normal smoothness, and object-mask opacity regularization. Stage II assigns albedo and roughness to each Gaussian, rasterizes albedo and roughness maps, uses rendered depth and normals to define per-pixel surface position and normal, samples incident directions, traces them to obtain $$\mathcal{G}_{i}(\mathbf{x}) = \exp\left(-\frac{1}{2}(\mathbf{x}-\mu_i)^\top \Sigma_i^{-1}(\mathbf{x}-\mu_i)\right), \qquad \Sigma_i = R_i S_i S_i^\top R_i^\top,$$1 and $$\mathcal{G}_{i}(\mathbf{x}) = \exp\left(-\frac{1}{2}(\mathbf{x}-\mu_i)^\top \Sigma_i^{-1}(\mathbf{x}-\mu_i)\right), \qquad \Sigma_i = R_i S_i S_i^\top R_i^\top,$$2, evaluates the rendering equation, and optimizes the final PBR color with regularizers (Gu et al., 2024).

The rendering equation is evaluated by Monte Carlo integration. For dielectric objects, IRGS uses stratified hemisphere sampling with $$\mathcal{G}_{i}(\mathbf{x}) = \exp\left(-\frac{1}{2}(\mathbf{x}-\mu_i)^\top \Sigma_i^{-1}(\mathbf{x}-\mu_i)\right), \qquad \Sigma_i = R_i S_i S_i^\top R_i^\top,$$3 incident directions per pixel, and estimates the PBR color as

$$\mathcal{G}_{i}(\mathbf{x}) = \exp\left(-\frac{1}{2}(\mathbf{x}-\mu_i)^\top \Sigma_i^{-1}(\mathbf{x}-\mu_i)\right), \qquad \Sigma_i = R_i S_i S_i^\top R_i^\top,$$4

To control cost, it imposes a per-iteration ray budget $$\mathcal{G}_{i}(\mathbf{x}) = \exp\left(-\frac{1}{2}(\mathbf{x}-\mu_i)^\top \Sigma_i^{-1}(\mathbf{x}-\mu_i)\right), \qquad \Sigma_i = R_i S_i S_i^\top R_i^\top,$$5 and reduces the number of evaluated pixels to $$\mathcal{G}_{i}(\mathbf{x}) = \exp\left(-\frac{1}{2}(\mathbf{x}-\mu_i)^\top \Sigma_i^{-1}(\mathbf{x}-\mu_i)\right), \qquad \Sigma_i = R_i S_i S_i^\top R_i^\top,$$6 (Gu et al., 2024). During relighting, it does not recursively trace full indirect bounces. Instead, it ray traces and alpha blends albedo, roughness, and normal to obtain aggregated values, then applies a split-sum approximation with a prefiltered environment map and sets $$\mathcal{G}_{i}(\mathbf{x}) = \exp\left(-\frac{1}{2}(\mathbf{x}-\mu_i)^\top \Sigma_i^{-1}(\mathbf{x}-\mu_i)\right), \qquad \Sigma_i = R_i S_i S_i^\top R_i^\top,$$7 (Gu et al., 2024).

6. Efficiency, empirical findings, and open technical questions

LiDAR-RT is implemented in PyTorch with custom CUDA kernels and uses NVIDIA OptiX for hardware-accelerated ray tracing. The reported training schedule is 30,000 iterations for Gaussian optimization and 500 epochs for the U-Net refinement, on a single RTX 4090 GPU, with near plane $$\mathcal{G}_{i}(\mathbf{x}) = \exp\left(-\frac{1}{2}(\mathbf{x}-\mu_i)^\top \Sigma_i^{-1}(\mathbf{x}-\mu_i)\right), \qquad \Sigma_i = R_i S_i S_i^\top R_i^\top,$$8 and ray-tracing chunk size 16. The reported throughput is about 20.1 FPS on Waymo and about 42.7 FPS on KITTI-360. The paper attributes efficiency to the BVH over compact proxy geometry, chunked sorted intersection processing, hardware ray tracing, front-to-back differentiable accumulation, and the avoidance of full global order storage in the backward pass (Zhou et al., 2024).

Its ablations clarify why the 2D formulation is used. Compared with 3D Gaussians, 2D Gaussians are reported to be slightly better in rendering quality and efficiency and to improve geometric metrics significantly. Additional ablations show that removing separate hit/drop modeling hurts quality, removing ray-drop refinement worsens results, removing ray-direction inputs to refinement loses dynamic-object detail, removing Chamfer loss degrades geometry, and removing normal-based initialization degrades scene reconstruction (Zhou et al., 2024). This suggests that, in LiDAR simulation, differentiable 2D Gaussian ray tracing is effective only when paired with explicit sensor-response modeling and geometric regularization.

IRGS is implemented in OptiX through PyTorch CUDA extensions. Because geometry changes during optimization, the BVH is updated every training iteration; the reported BVH update cost is about 3 ms. The implementation uses $$\mathcal{G}_{i}(\mathbf{x}) = \exp\left(-\frac{1}{2}(\mathbf{x}-\mu_i)^\top \Sigma_i^{-1}(\mathbf{x}-\mu_i)\right), \qquad \Sigma_i = R_i S_i S_i^\top R_i^\top,$$9 for per-ray sorting, stops ray tracing when transmittance drops below 0.03, uses $\mu_i$0 rays for rendering-equation evaluation, sets the maximum $\mu_i$1 per iteration, renders a complete image at $\mu_i$2 in about 1 second, and trains in about 40 minutes on one RTX 3090 (Gu et al., 2024).

IRGS ablations isolate the role of differentiable tracing of indirect light. Detaching the indirect term degrades results, especially for indirect radiance and albedo estimation. Removing indirect during training harms albedo estimation, and removing indirect during relighting reduces relighting quality. The number of rays is also consequential: $\mu_i$3 is poor, $\mu_i$4 is better, and $\mu_i$5 is best. The reported qualitative effects are more plausible indirect illumination, cleaner occlusion, more realistic relit images, and better material separation (Gu et al., 2024).

A recurrent technical question in this literature is whether Gaussian ray tracing should operate on volumetric 3D kernels or on 2D surface-like disks. VoGE demonstrates that differentiable ray tracing over 3D Gaussian ellipsoids can be efficient and useful for analysis-by-synthesis, but it also relies on approximate closed-form aggregation, truncation to the nearest $\mu_i$6 kernels, and assumptions such as a perspective camera (Wang et al., 2022). IRGS argues that 3D Gaussian ray tracing is problematic for inverse rendering because the ray–splat intersection is ambiguous and tracing a pretrained 3DGS checkpoint causes noticeable rendering degradation. LiDAR-RT reports that 2D Gaussians are slightly better than 3D Gaussians in both rendering quality and efficiency for LiDAR re-simulation (Gu et al., 2024). A plausible implication is that surface-like 2D Gaussian primitives are particularly advantageous when the target task requires a geometrically meaningful ray hit, such as active sensing or the evaluation of visibility and indirect radiance in the rendering equation.