3D Gaussian Splatting

Lesson 09's NeRF is beautiful and slow: an implicit radiance field queried by a neural network hundreds of times per ray, so a single frame takes seconds. This lesson keeps the one idea that made NeRF trainable — differentiable front-to-back compositing — but swaps the representation out from under it. Instead of an implicit field you march, we place explicit primitives you rasterize. That single swap buys real-time rendering (100+ FPS) and easy editing, and in doing so it resolves the explicit-vs-implicit tension from lesson 01 decisively in favor of speed.

1 · Why NeRF is slow, and what we want instead

Recall the NeRF forward pass from lesson 09. To color one pixel you cast a ray, pick tens to low-hundreds of sample points along it, and at every sample evaluate an MLP F_θ(x, dir) → (color, density) to get the emitted color and volume density; then you integrate those samples front-to-back. Multiply out the arithmetic: (pixels per frame) × (samples per ray) × (one full MLP forward per sample). For a 800×800 image at ~192 samples/ray that is on the order of 10⁸ network evaluations per frame. Even with the hash-grid acceleration of Instant-NGP, the representation is implicit — there is no geometry to draw, only a field to interrogate sample-by-sample — so rendering is fundamentally ray marching through a network. That is why vanilla NeRF renders in seconds per frame and trains in hours.

Now write down what we actually want for a product that has to render and be edited interactively:

Each term is doing work. Explicit primitives means the scene is a list of things with coordinates — no per-sample network to evaluate. Rasterization means we project those primitives to the screen and let the GPU's native strength (drawing lots of small overlapping shapes) do the rendering, instead of marching rays. Differentiability means we keep the property from lesson 01 §4 that made the whole inverse-rendering loop work: gradients flow from the image error back into the primitives, so we can still fit the scene to photos by gradient descent.

2 · The primitive — a 3D Gaussian

A 3DGS scene is a set of a few million 3D Gaussians, each a small anisotropic blob. One Gaussian carries four things:

• Mean μ ∈ ℝ³ — where the blob sits in world space (its position).
• Covariance Σ ∈ ℝ³ˣ³ — the blob's size, its anisotropic stretch, and its orientation, all in one matrix. This is the interesting parameter, and §2 is mostly about how to store it safely.
• Opacity α ∈ [0,1] — how much the blob occludes what is behind it (its peak alpha).
• View-dependent color via spherical-harmonic (SH) coefficients — instead of one RGB, each Gaussian stores a small set of SH coefficients per channel, and the color emitted toward the camera is the SH series evaluated in the viewing direction. Degree-0 SH is a flat color; adding higher-degree bands lets the blob look different from different angles, which is how 3DGS reproduces specular highlights and glossy reflections — the same view-dependence NeRF got from feeding dir into its MLP.

The blob's spatial influence at a point x is an unnormalized Gaussian centered at μ with covariance Σ:

The quadratic form (x-μ)ᵀ Σ^-1 (x-μ) is the squared Mahalanobis distance — the "distance in units of the blob's own spread." It equals 1 on the 1σ ellipsoid, so Σ literally is the shape: its eigenvectors are the ellipsoid's axes and the square roots of its eigenvalues are the axis half-lengths.

The covariance must stay valid — enter quaternions (lesson 02)

Here is the subtlety that trips people up. A covariance matrix must be symmetric positive semidefinite (PSD) — otherwise Σ^-1 and the quadratic form are meaningless (a negative eigenvalue would make G(x) blow up along some direction instead of decaying). A symmetric 3×3 has 6 independent entries, and it is tempting to just make those 6 numbers learnable and let the optimizer move them. Do not.

The construction factors the covariance into a scale and a rotation. Let S = \mathrm{diag}(s_x, s_y, s_z) be a diagonal scale (the blob's extent along its own three axes) and let R be a rotation matrix built from a unit quaternion q. Define

Why this is guaranteed valid: for any vector v, vᵀΣv = vᵀ R S Sᵀ Rᵀ v = ‖Sᵀ Rᵀ v‖² ≥ 0, so Σ is PSD automatically — no constraint to enforce, no re-projection step. The factorization also cleanly separates the two things you want to control independently: S is the shape/size (how stretched the blob is along its own axes) and R is the orientation (which way those axes point). You optimize the three log-scales \log s (log-space keeps them positive) and the four quaternion components (renormalized to unit length each step), and reconstruct Σ from them for rendering.

3 · Rendering by splatting — the EWA covariance transform

To draw the scene we must get each 3D Gaussian onto the 2D image. "Splatting" means projecting a 3D Gaussian to a 2D Gaussian in image space and stamping (splatting) that 2D footprint onto the pixels. Two pieces project: the mean and the covariance.

The mean projects like any 3D point — through the world→camera rigid transform and then the pinhole projection of CV lesson 04: μ goes to camera coordinates, then to a pixel via the intrinsics K. Nothing new.

The covariance is the interesting part. A 2D Gaussian in image space also has a 2×2 covariance Σ', and it is obtained by pushing the 3D covariance through the same transformation the mean underwent — but a covariance transforms by sandwiching with the transformation's linear part. The EWA (Elliptical Weighted Average) splatting result is:

where W is the world→camera rotation (the viewing transform's linear part — it reorients the covariance into camera coordinates) and J is the Jacobian of the projective map evaluated at μ. The sandwich form is the standard rule for how a covariance transforms under a linear map y = My: \mathrm{Cov}(y) = M\,\mathrm{Cov}(x)\,Mᵀ. The composition J W is that linear map: first rotate into the camera frame (W), then apply the local linearization of the perspective projection (J).

The upshot: after splatting, each visible Gaussian is a 2D elliptical bump on the image with center = projected μ and shape = Σ'. Evaluating that 2D Gaussian at a pixel gives a scalar in [0,1] — the blob's spatial weight at that pixel — which we combine with the blob's opacity next.

4 · Differentiable tile-based rasterization

Now blend the splats into a picture, differentiably and fast. The 3DGS rasterizer does it in three stages:

1. Tile. Partition the image into small tiles (e.g. 16×16 pixels). For each Gaussian, find which tiles its 2D footprint Σ' overlaps, and add it to those tiles' lists. This turns "which of the millions of Gaussians touch this pixel?" into a cheap local lookup and lets each tile be rasterized independently (great for the GPU).
2. Sort by depth. Within the visible set, sort the Gaussians front-to-back by their camera-space depth. This is one global per-view sort (a fast GPU radix sort keyed by depth + tile), not a per-pixel sort — every pixel in a tile then walks the same depth-ordered list. Correct compositing requires a consistent front-to-back order, and this sort supplies it.
3. Blend front-to-back, per pixel. Walk the sorted list and accumulate color with alpha compositing, stopping early once the pixel is opaque (transmittance ≈ 0).

The blend is the alpha-over compositing equation:

where c_i is the Gaussian's SH color for this view, α_i is its learned opacity, and G'_i(\text{pixel}) is its projected 2D Gaussian evaluated at that pixel. The product \prod_{j<i}(1-α'_j) is the accumulated transmittance T_i — the fraction of light that survived all the nearer, partially-opaque blobs — so each blob contributes its color scaled by how visible it still is.

Why it is differentiable. Every operation here — the Gaussian G'_i, the product α_i\cdot G'_i, the running transmittance, the weighted sum — is smooth in the Gaussians' parameters. There is no hard coverage test (the discontinuity that lesson 01 §4 warned breaks mesh rasterization); a Gaussian's soft falloff gives every pixel a nonzero gradient with respect to the blob's mean, covariance, opacity, and color. So the rasterizer is a differentiable renderer R in exactly the sense the inverse-rendering loop needs, and it runs on the GPU at interactive rates.

5 · Optimization and adaptive density control

Fitting a scene is the loop from lesson 01, instantiated. Given a set of posed training images, render each view with the rasterizer above and compare to the ground-truth photo with a photometric loss — the same family NeRF used: an L1 term plus a D-SSIM (structural-similarity) term that penalizes local structural error, not just per-pixel color:

Because the whole forward pass is differentiable, gradients flow to every parameter of every Gaussian: the mean μ, the log-scales, the orientation quaternion, the opacity α, and the SH color coefficients. A plain optimizer (Adam) descends all of them jointly.

But how many Gaussians, and where? You do not know the right count up front, so 3DGS changes the set itself during training — this is adaptive density control, and it is what lets a few million blobs organize into a sharp scene:

• Initialize from the SfM point cloud. The camera poses come from Structure-from-Motion (CV lesson 05), and SfM also produces a sparse 3D point cloud as a by-product. 3DGS seeds one small Gaussian at each such point — a far better start than random, because those points already sit on real surfaces. (This is the same SfM initialization NeRF-style pipelines rely on, reused here.)
• Clone in under-reconstructed regions. Where a region is missing detail, the tell-tale is a small Gaussian carrying a large positional gradient (the optimizer wants to move it to cover more than it can). 3DGS clones such Gaussians — duplicating the blob and nudging the copy along the gradient — to add coverage where the scene is under-populated.
• Split in over-reconstructed regions. Where a single Gaussian has grown too large and is smearing detail (large scale, large positional gradient), 3DGS splits it into two smaller Gaussians (scaled down and offset), trading one coarse blob for finer ones that can capture high-frequency structure.
• Prune the useless. Gaussians whose opacity has decayed toward zero (α → 0) contribute nothing, so they are periodically removed, keeping the count — and memory and render cost — in check. (Opacity is also reset low every so often to force genuinely-needed blobs to re-earn their opacity, which lets floaters be pruned.)

The result is an automatic allocation of capacity: dense clusters of tiny Gaussians on detailed, textured surfaces; a few big ones on flat walls; nothing wasted in empty space. This grow/prune behavior is the explicit-representation counterpart of NeRF letting an MLP implicitly distribute its capacity — here the capacity is a literal, editable list of primitives you can add to and delete from.

6 · NeRF vs 3DGS — the comparison

Both fit a scene to posed images by differentiable front-to-back compositing under a photometric loss. Everything else differs, and the differences all trace back to the one choice — implicit field vs explicit primitives:

The widget below is 3DGS in 2D: a handful of 2D Gaussians, each with a covariance you build from a rotation angle and two scales — the planar version of §2's Σ = RSSᵀRᵀ — rasterized by the exact front-to-back blend of §4. Stretch and rotate one blob to cover a diagonal streak with a single anisotropic Gaussian; reorder the depths and watch which blob wins the pixel ("popping"); dial opacity to feel the difference between hard occlusion and see-through blending.

// Per-Gaussian covariance from rotation θ and scales (s_x, s_y):  Σ = R diag(s_x², s_y²) Rᵀ
// We precompute Σ⁻¹ (a 2×2 inverse) once per Gaussian, not per pixel.
function invCov(sx, sy, th){
  var c = Math.cos(th), s = Math.sin(th);
  // Σ = [[a,b],[b,d]]
  var a = c*c*sx*sx + s*s*sy*sy;
  var d = s*s*sx*sx + c*c*sy*sy;
  var b = c*s*(sx*sx - sy*sy);
  var det = a*d - b*b;
  return { ia: d/det, ib: -b/det, id: a/det };   // Σ⁻¹
}

// RASTERIZE: iterate pixels (coarse step for speed); blend blobs FRONT-TO-BACK by depth.
for (var py = 0; py < H; py += STEP)
for (var px = 0; px < W; px += STEP){
  var T = 1, r = bgR, g = bgG, b = bgB;            // transmittance + accumulated colour
  for (var k = 0; k < order.length; k++){          // order = indices sorted by depth (near→far)
    var G = gauss[order[k]], dx = px - G.mx, dy = py - G.my;
    var m2 = G.ia*dx*dx + 2*G.ib*dx*dy + G.id*dy*dy;   // Mahalanobis²
    var val = Math.exp(-0.5*m2);                        // 2D Gaussian weight ∈ (0,1]
    var aP = G.opacity * val;                           // α′ = opacity · g
    r += T*aP*G.r; g += T*aP*G.g; b += T*aP*G.b;        // C += T·α′·colour
    T *= (1 - aP);                                       // update transmittance
    if (T < 0.003) break;                                // pixel opaque → early out
  }
  fillRect(px, py, STEP, STEP, r, g, b);
}
// Reordering depths changes `order` → a different blob blends on top ("popping").

Where this points next

Both branches of neural rendering — implicit NeRF and explicit 3DGS — now share a spine: fit one scene to its own posed photos by differentiable compositing. But notice what they never do: generalize. Optimize a NeRF or a splat cloud for a room, and you learn that room and nothing else; a new scene means training again from scratch. That is fine for capture, useless for a model that should look at a novel object and infer its 3D shape. Lesson 11 makes the turn to learning across scenes: neural networks that operate directly on 3D data — point clouds (PointNet and permutation invariance), voxels and sparse convolution (the O(N³) escape from lesson 01 §3), and meshes — so a network trained on many objects can predict on ones it has never seen. The representations from lesson 01 return, this time as the inputs and outputs of learned models rather than things you optimize one at a time.

Interview prompts

• Why is 3DGS real-time when NeRF is not, given they share the same compositing equation? (§1, §4 — NeRF ray-marches an implicit field with hundreds of MLP queries per ray; 3DGS rasterizes explicit primitives and blends a handful of depth-sorted splats per pixel — no per-sample network, hardware-friendly sorted rasterization.)
• Why store a Gaussian's covariance as Σ = RSSᵀRᵀ with a quaternion, instead of its 6 raw entries? (§2 — free entries drift to a negative eigenvalue → indefinite, invalid Σ; the factorization is PSD by construction and separates scale from orientation, reusing lesson 02's quaternion.)
• Derive how a 3D covariance projects to 2D, and say why it is only approximate. (§3 — Σ' = JWΣWᵀJᵀ by the covariance-sandwich rule; J is the local affine linearization of the nonlinear perspective map at μ, so large / far-off-axis Gaussians are distorted.)
• Write the front-to-back blending equation and relate it to NeRF. (§4 — C = \sum_i c_iα'_i\prod_{j<i}(1-α'_j) with α'_i = α_iG'_i; identical alpha compositing to NeRF §3, summed over sorted splats not ray samples.)
• What does adaptive density control do, and how is the optimization initialized? (§5 — clone under-reconstructed (small blob, large positional grad), split over-large blobs, prune α→0; initialized from the SfM sparse point cloud of CV lesson 05.)
• When would you choose 3DGS over NeRF, and what artifacts does each show? (§6 — 3DGS for real-time rendering/editing if memory allows (hundreds of MB, popping + spiky blobs); NeRF for compactness / a continuous field (small weights, floaters).)