3D Gaussian Splatting for Real-Time Radiance Field Rendering

The tension the paper resolves

By 2023 novel-view synthesis had split into two camps, each forced to give something up:

• Continuous / neural representations (NeRF and descendants). Optimize an MLP queried by volumetric ray-marching. Continuity aids optimization and gives beautiful results, but rendering requires stochastic sampling of many points per ray — slow and noisy. The quality leader, Mip-NeRF360, needed up to 48 hours of training and rendered at seconds per frame.
• Fast grid / hash methods (InstantNGP, Plenoxels). Storing interpolatable features in voxel/hash grids cut training to minutes and reached interactive (~10–15 fps) rates, but they still ray-march, still struggle with empty space, and their quality is capped by the structured acceleration grid.

Nobody could get all three at once — SOTA quality, fast training, and real-time 1080p rendering of full scenes. That is the gap 3DGS closes.

The core insight

Meshes and points map cleanly onto fast GPU rasterization. NeRFs win on quality because volumetric α-blending is a forgiving, gradient-friendly image-formation model. 3DGS asks: what primitive is simultaneously (a) a differentiable volumetric element and (b) trivially and cheaply rasterizable?

The answer is an anisotropic 3D Gaussian — a smooth differentiable volumetric bump (so optimization behaves like a radiance field) that projects to a 2D Gaussian “splat” blendable with ordinary α-compositing (so rendering behaves like point rasterization). Empty space contains no Gaussians, so no compute is wasted there.

The three contributions

1. Anisotropic 3D Gaussians as a high-quality, unstructured representation of radiance fields.
2. An optimization method for Gaussian properties, interleaved with adaptive density control.
3. A fast, differentiable, visibility-aware GPU renderer supporting anisotropic splatting and efficient backpropagation.

Input is the same as NeRF’s: SfM-calibrated cameras, which produce a sparse point cloud for free. Unlike most prior point-based methods, 3DGS needs only SfM points — no Multi-View Stereo. On synthetic data it even works from random initialization. Converged scenes hold 1–5 million Gaussians.

Traditional reconstruction and rendering

Early novel-view synthesis used light fields. Structure-from-Motion (SfM) enabled turning photo collections into novel views, producing a sparse point cloud as a by-product. Multi-View Stereo (MVS) gave full 3D reconstruction feeding image-based rendering methods that re-project and blend inputs — but these fail on unreconstructed regions and MVS “over-reconstruction” (invented geometry).

Neural rendering and radiance fields

NeRF introduced positional encoding and importance sampling with a large MLP — great quality, slow speed. The quality SOTA became Mip-NeRF360. Acceleration methods exploit three levers: spatial data structures for features, cheaper encodings, and smaller/absent MLPs:

• InstantNGP — multiresolution hash grid + occupancy grid + small MLP.
• Plenoxels — sparse voxel grid interpolating a continuous density field, no neural network at all.

Both use spherical harmonics for directional appearance; both struggle with empty space and are limited by structured grids.

Point-based rendering and the shared image-formation model

Naïve point rendering has holes and aliasing, fixed by “splatting” primitives larger than a pixel (the EWA lineage of Zwicker et al.). The key theoretical point: point-based α-blending and NeRF volumetric rendering share the same image-formation model.

Volumetric rendering along a ray:

\[ C = \sum_{i=1}^{N} T_i\,\bigl(1 - \exp(-\sigma_i \delta_i)\bigr)\,\mathbf{c}_i, \qquad T_i = \exp\!\left(-\sum_{j=1}^{i-1}\sigma_j \delta_j\right). \tag{1}\]

Re-writing with \(\alpha_i = 1 - \exp(-\sigma_i \delta_i)\):

\[ C = \sum_{i=1}^{N} T_i\,\alpha_i\,\mathbf{c}_i, \qquad T_i = \prod_{j=1}^{i-1}(1-\alpha_j). \tag{2}\]

A typical neural point method blends \(N\) ordered points overlapping a pixel:

\[ C = \sum_{i \in N} \mathbf{c}_i\,\alpha_i \prod_{j=1}^{i-1}(1-\alpha_j), \tag{3}\]

where \(\alpha_i\) evaluates a 2D Gaussian with covariance \(\Sigma\) times a learned per-point opacity. Equation 2 and Equation 3 are the same model; the difference is algorithmic — NeRF implicitly represents space and pays for random sampling, while points are explicit, unstructured, and editable.

Pulsar inspired the tile-based sorted renderer, but is order-independent. 3DGS keeps approximate conventional α-blending on sorted splats (visibility respected), back-propagates to all splats in a pixel, and rasterizes anisotropic splats — the ingredients behind its quality. Prior 3D-Gaussian uses targeted single isolated objects; 3DGS scales to complete, complex, high-depth-complexity scenes, indoor and outdoor.

Why Gaussians, and how they are defined

The goal: a primitive that inherits differentiable-volumetric properties yet is unstructured and explicit for fast rendering — and that works from normal-free SfM points. 3D Gaussians fit. A Gaussian is defined by a full 3D covariance \(\Sigma\) in world space, centered at mean \(\mu\):

\[ G(x) = e^{-\frac{1}{2}(x-\mu)^{T}\,\Sigma^{-1}\,(x-\mu)}. \tag{4}\]

In blending, \(G(x)\) is multiplied by the Gaussian’s opacity \(\alpha\).

Projecting 3D Gaussians to 2D (EWA splatting)

For rendering, each Gaussian becomes a 2D splat. Following Zwicker et al.’s EWA framework, given a viewing transform \(W\), the camera-space covariance is:

\[ \Sigma' = J\,W\,\Sigma\,W^{T}\,J^{T}, \tag{5}\]

where \(J\) is the Jacobian of the affine approximation of the projective transform. Dropping the third row and column of \(\Sigma'\) yields a \(2\times2\) covariance — a proper 2D Gaussian splat.

A parameterization that stays valid under gradient descent

You cannot optimize the raw entries of \(\Sigma\) directly: a covariance is only physically meaningful when positive semi-definite, and unconstrained gradient steps produce invalid matrices. Since \(\Sigma\) describes an ellipsoid, factor it into a rotation and a scale:

\[ \Sigma = R\,S\,S^{T}\,R^{T}. \tag{6}\]

Stored separately and optimized independently: a 3D scale vector \(s\) (→ diagonal \(S\)) and a quaternion \(q\) (normalized → rotation \(R\)). As expressive as a raw covariance, but always valid. To avoid autodiff overhead on this hot path, all gradients are derived explicitly (Section 1.12).

The full per-Gaussian parameter set

The optimization step

Iterate: render a training view, compare to the captured image, back-propagate, step. Because 3D→2D projection is ambiguous, the optimizer must create, destroy, and move geometry. Key choices:

• Optimizer: Adam SGD, with custom CUDA kernels on the critical path (rasterization is the bottleneck).
• Activations: sigmoid on \(\alpha\) (keeps it in \([0,1)\) with smooth gradients); exponential on covariance scale.
• Initial covariance: isotropic, axis length = mean distance to the 3 nearest SfM points.
• LR schedule: exponential decay (à la Plenoxels), applied to positions only.
• Loss: \(\mathcal{L}_1\) + D-SSIM:

\[ \mathcal{L} = (1-\lambda)\,\mathcal{L}_1 + \lambda\,\mathcal{L}_{\text{D-SSIM}}, \qquad \lambda = 0.2. \tag{7}\]

Adaptive control of Gaussians — clone, split, prune

Every 100 iterations (after warm-up), density is adjusted. Two failure modes are targeted, and both show up as large view-space positional gradients (the optimizer straining to move Gaussians):

• Under-reconstruction — missing geometry; small Gaussians can’t cover a region.
• Over-reconstruction — one big Gaussian sprawling over detailed geometry.

Gaussians with average view-space position-gradient magnitude above \(\tau_{pos} = 0.0002\) are densified:

• Under-reconstruction → CLONE: copy the small Gaussian, shift the copy along the positional gradient. Adds volume and a Gaussian.
• Over-reconstruction → SPLIT: replace the large Gaussian with two, dividing scale by \(\phi = 1.6\), sampling positions using the original as a PDF. Conserves volume, adds a Gaussian.

Pruning / count control. Remove Gaussians with \(\alpha < \epsilon_\alpha\). To fight near-camera “floaters,” every N = 3000 iterations \(\alpha\) is reset near zero; the optimizer re-raises \(\alpha\) only where needed, and culling removes the rest. Oversized world-space or viewspace Gaussians are periodically removed. Gaussians stay Euclidean primitives at all times — no space warping/contraction (unlike Mip-NeRF360 or Plenoxels).

Densification logic (reconstruction of Fig. 4 / Algorithm 1):

flowchart TD
    A["Refinement iteration<br/>(every 100 iters)"] --> B{"For each Gaussian"}
    B --> C{"alpha &lt; eps_a<br/>or too large?"}
    C -- "yes" --> D["Prune / remove"]
    C -- "no" --> E{"View-space pos.<br/>gradient &gt; tau_pos?"}
    E -- "no" --> F["Leave unchanged"]
    E -- "yes" --> G{"Scale ||S|| &gt; tau_S?"}
    G -- "large<br/>(over-reconstruction)" --> H["SPLIT into 2<br/>scale / phi (1.6)<br/>sample new positions"]
    G -- "small<br/>(under-reconstruction)" --> I["CLONE<br/>copy + move along<br/>positional gradient"]

The overall training loop (Algorithm 1):

flowchart TD
    A["M = SfM points (positions)"] --> B["S, C, A = InitAttributes()<br/>covariances, colors, opacities"]
    B --> C["i = 0"]
    C --> D{"not converged?"}
    D -- "no" --> Z["Done then trained model"]
    D -- "yes" --> E["V, GT = SampleTrainingView()"]
    E --> F["I = Rasterize(M,S,C,A,V)  (Alg. 2)"]
    F --> G["L = Loss(I, GT)"]
    G --> H["M,S,C,A = Adam(grad L)"]
    H --> I{"IsRefinementIteration(i)?"}
    I -- "yes" --> J["Prune (alpha&lt;eps or too large)<br/>+ Densify (split / clone)"]
    I -- "no" --> K["i = i + 1"]
    J --> K
    K --> D

Tile-based forward pass

1. Tile the screen into 16×16 pixel tiles.
2. Cull against the view frustum — keep Gaussians whose 99% confidence ellipsoid intersects it; a guard band trivially rejects unstable extreme-position Gaussians.
3. Instantiate each Gaussian once per overlapped tile, with a 64-bit key: low 32 bits = view-space depth, high bits = tile ID.
4. Sort all instances with a single GPU radix sort. Whole-image depth order resolved at once — no per-pixel sorting. Blending is therefore approximate, but error vanishes as splats approach pixel size.
5. Build per-tile ranges (one thread per key, comparing high bits with neighbors).
6. Rasterize: one thread block per tile; threads collaboratively load Gaussians into shared memory; each pixel walks its list front-to-back, accumulating color and \(\alpha\), stopping when \(\alpha\) saturates; the tile terminates when all pixels saturate.

Backward pass without storing per-pixel lists

• Reuse the sorted array and per-tile ranges, traversing back-to-front.
• A pixel processes a point only if its depth ≤ the depth of the last point that contributed in the forward pass.
• Instead of storing every intermediate opacity, each pixel stores only the total accumulated \(\alpha\), and the backward pass divides by each point’s \(\alpha\) to recover the coefficients.

Because the only stopping criterion is \(\alpha\)-saturation (not a fixed count), the method handles arbitrary depth complexity with no scene-specific tuning — critical to quality (the ablation capping gradient depth to 10 splats lost ~11 dB PSNR).

Rasterizer (Algorithm 2):

flowchart TD
    A["RASTERIZE(w, h, M, S, C, A, V)"] --> B["CullGaussian(p, V)<br/>frustum + guard band"]
    B --> C["M', S' = ScreenspaceGaussians(M,S,V)<br/>EWA projection to 2D"]
    C --> D["T = CreateTiles(w, h)  (16x16)"]
    D --> E["L, K = DuplicateWithKeys(M', T)<br/>key = tileID and depth"]
    E --> F["SortByKeys(K, L)<br/>single global GPU radix sort"]
    F --> G["R = IdentifyTileRanges(T, K)"]
    G --> H["for each tile t, each pixel i:<br/>BlendInOrder front-to-back<br/>until alpha saturates"]
    H --> I["return image I"]

Numerical stability

Reconstructing opacities by division risks divide-by-zero. Safeguards: skip blending with \(\alpha < \epsilon\) (\(\epsilon = \tfrac{1}{255}\)), clamp \(\alpha \le 0.99\), and — before adding a Gaussian — stop front-to-back blending if accumulated opacity would exceed 0.9999.

Implementation details

• Stack: Python + PyTorch; custom CUDA kernels for rasterization; NVIDIA CUB radix sort; interactive viewer on SIBR.
• Resolution warm-up: start at ¼ resolution, upsample after 250 and 500 iterations.
• Progressive SH: optimize only the zero-order (diffuse) band first, add one band every 1000 iterations up to 4 bands — prevents diffuse-color corruption on partial-hemisphere captures.
• Hardware: single NVIDIA A6000 (Mip-NeRF360 trained on A100s).

Quantitative results — real scenes (Table 1)

13 real scenes across Mip-NeRF360, Tanks & Temples, Deep Blending; same hyperparameters; every 8th photo held out; SSIM ↑, PSNR ↑, LPIPS ↓.

Ours-30K matches or beats Mip-NeRF360’s quality while training in ~35–45 min (vs. 48 h) and rendering at 130–150+ fps (vs. 0.06–0.14 fps). Ours-7K (≈5–7 min) already surpasses InstantNGP/Plenoxels on most metrics at 160–197 fps. The cost is memory: hundreds of MB vs. tens of MB for grid/hash methods.

Quantitative results — synthetic Blender (Table 2)

From 100K random Gaussians (no SfM), pruned to 6–10K meaningful ones, final 200–500K, rendered at 180–300 fps.

3DGS posts the best average PSNR (33.32) even from random init.

Compactness. Against Zhang et al. 2022, 3DGS breaks even on PSNR in 2–4 minutes using ~¼ the point count, at ~3.8 MB vs. 9 MB (2 SH degrees for fairness).

Ablations (Table 3)

• Initialization from SfM. Random init still avoids total failure but degrades in the background and leaves un-removable floaters (Fig. 7).
• Densification: clone vs. split. Split helps background reconstruction; clone aids thin structures and convergence (Fig. 8).
• Unlimited gradient depth. Capping gradients to N = 10 front-most splats caused unstable optimization and −11 dB PSNR on Truck (Fig. 9).
• Anisotropic covariance. Replacing full covariance with a single isotropic radius significantly lowers quality (Fig. 10).
• Spherical harmonics. SH improves PSNR by modeling view-dependent effects.

Limitations

• Artifacts in poorly observed regions (elongated/“splotchy” Gaussians) — Mip-NeRF360 also struggles (Figs. 11, 12).
• Popping artifacts for large Gaussians, aggravated by guard-band rejection and simple global-sort visibility; anti-aliasing and principled culling are future work.
• No regularization yet.
• Memory: peak training memory can exceed 20 GB (unoptimized prototype); models are hundreds of MB + 30–500 MB rasterizer overhead.
• Very large scenes may need a lower position learning rate.