3D Gaussian Splatting: the rasterizer

3D Gaussian Splatting (Kerbl, Kopanas, Leimkühler & Drettakis, SIGGRAPH 2023) represents a scene as a few million fuzzy, coloured, semi-transparent ellipsoids and renders them by rasterization — no network, no ray marching — at 100+ fps at 1080p, with quality that matched or beat the best NeRFs of the time. The renderer is also differentiable, which is how the ellipsoids are fitted to photographs (next skill, gfx-differentiable-rendering).

This lesson walks through the forward pass of the reference CUDA rasterizer (graphdeco-inria/diff-gaussian-rasterization) line by line, and notes where gsplat (nerfstudio) does things differently. You'll need quaternions (gfx-quaternions), the pinhole camera (gfx-cameras-projection), front-to-back compositing (gfx-volume-rendering) and degree-0 SH colour (gfx-spherical-harmonics).

1. What a splat is

Each Gaussian carries:

Storing the unconstrained version of each quantity is a recurring trick: the optimiser can push a raw number anywhere and the activation keeps the physical value legal. In the reference those activations happen in Python (torch.exp, torch.sigmoid, F.normalize) before the CUDA kernel sees the data; in the build task your preprocess applies them.

2. The 3D covariance: Σ = R S Sᵀ Rᵀ

A Gaussian's shape is its covariance Σ: G(x) = exp(−½ (x − μ)ᵀ Σ⁻¹ (x − μ)). Σ must be symmetric positive semi-definite. If you optimised its six entries directly, a gradient step could make it indefinite — a negative variance. So 3DGS factors it as a scaled-then-rotated unit sphere:

M = R S, Σ = M Mᵀ = R S Sᵀ Rᵀ, S = diag(s), R = R(q/|q|).

For any s and q, xᵀ Σ x = |Mᵀ x|² ≥ 0: valid by construction. The order matters — scale first, then rotate. (S R)(S R)ᵀ = S R Rᵀ S = S², and the rotation vanishes.

#include <array>
#include <cmath>
#include <numbers>
#include <print>

using Mat3 = std::array<std::array<float, 3>, 3>;

// Rotation matrix of a unit quaternion (w, x, y, z).
Mat3 rotation(float w, float x, float y, float z) {
    const float n = std::sqrt(w * w + x * x + y * y + z * z);
    w /= n; x /= n; y /= n; z /= n;
    return {{{1 - 2 * (y * y + z * z), 2 * (x * y - w * z), 2 * (x * z + w * y)},
             {2 * (x * y + w * z), 1 - 2 * (x * x + z * z), 2 * (y * z - w * x)},
             {2 * (x * z - w * y), 2 * (y * z + w * x), 1 - 2 * (x * x + y * y)}}};
}

// Σ_ij = Σ_k R_ik s_k² R_jk  — i.e. R S Sᵀ Rᵀ, without forming S.
Mat3 cov3d(std::array<float, 3> s, const Mat3& R) {
    Mat3 out{};
    for (int i = 0; i < 3; ++i)
        for (int j = 0; j < 3; ++j)
            for (int k = 0; k < 3; ++k) out[i][j] += R[i][k] * s[k] * s[k] * R[j][k];
    return out;
}

int main() {
    const float half = std::numbers::pi_v<float> / 8;   // a 45° turn about z
    const Mat3 R = rotation(std::cos(half), 0, 0, std::sin(half));
    const Mat3 S = cov3d({2.0f, 0.5f, 1.0f}, R);
    for (const auto& row : S) std::println("{:8.4f}{:8.4f}{:8.4f}", row[0], row[1], row[2]);
}

A 2 × 0.5 ellipse turned 45° gets equal x and y variances (4 + 0.25)/2 and a strong correlation (4 − 0.25)/2 off the diagonal.

Reading the reference. computeCov3D writes glm::mat3 M = S * R; Sigma = transpose(M) * M; — which looks like the S R bug. It isn't: glm is column-major, and the glm::mat3(…) constructor fills columns, so the matrix the code builds row by row is Rᵀ. Then M = S Rᵀ and Mᵀ M = R S S Rᵀ. When you port CUDA/glm code, check every transpose against the math, not against the spelling.

3. EWA projection to the screen

Zwicker et al.'s EWA splatting (2001) observed that an affine transform of a Gaussian is again a Gaussian: x ↦ A x + b maps Σ to A Σ Aᵀ. Perspective isn't affine, so it is linearised at the Gaussian's mean. With the view rotation W and camera-space mean t = W μ + t_cam, the projection (x, y, z) ↦ (fx·x/z, fy·y/z) has Jacobian

J = [[fx/z, 0, −fx·x/z²], [0, fy/z, −fy·y/z²], [0, 0, 0]],

and the screen-space covariance is the upper-left 2×2 of Σ′ = J W Σ Wᵀ Jᵀ.

Two details from computeCov2D that matter in practice:

• The clamp. Before building J, the reference clamps t.x/t.z and t.y/t.z to ±1.3·tan(fov/2), then multiplies back by t.z. The linearisation is poor far off-axis, and a Gaussian just outside the view would otherwise get a huge, smeared footprint across the screen. The 2D mean is not clamped.
• The dilation. cov[0][0] += 0.3f; cov[1][1] += 0.3f; — convolving two Gaussians adds their covariances, so this is a screen-space low-pass filter of σ ≈ 0.55 px: every splat becomes at least about a pixel wide, sub-pixel splats don't alias, and Σ′ stays invertible. The conic, the radius, everything downstream uses the dilated Σ′.

Aside: Mip-Splatting. The dilation has a cost: zoom out and every splat is fattened by the same 0.3 px², so small objects bloat and thin structures brighten. Mip-Splatting (Yu et al., CVPR 2024) adds a 3D smoothing filter bound to the training views' sampling rate and replaces the dilation with a "2D Mip filter" that scales opacity by √(det Σ′ / det(Σ′ + 0.3 I)) so the filter conserves energy. The reference rasterizer later gained exactly this as its antialiasing flag, and gsplat as rasterize_mode="antialiased".

The renderer evaluates the Gaussian with the inverse, the conic [[a b] [b c]]⁻¹ = [[c −b] [−b a]] / (ac − b²), and needs a bound on its extent: 3σ along the major axis. The eigenvalues of a symmetric 2×2 are λ = mid ± √(mid² − det) with mid = (a + c)/2; the reference floors the term under the root at 0.1 and takes radius = ceil(3·√λ_max) pixels.

#include <algorithm>
#include <cmath>
#include <print>

int main() {
    // A camera-space mean and an isotropic 3D covariance (σ = 0.1), camera at the origin.
    const float fx = 100, fy = 100, x = 0.4f, y = -0.2f, z = 2.0f, var = 0.01f;
    // J at the mean (W = I here, and x/z = 0.2 is well inside the 1.3·tan(fov/2) clamp).
    const float j00 = fx / z, j02 = -fx * x / (z * z), j11 = fy / z, j12 = -fy * y / (z * z);
    // Σ' = J Σ Jᵀ with Σ = var·I, upper-left 2×2:
    float a = var * (j00 * j00 + j02 * j02);
    float b = var * (j02 * j12);
    float c = var * (j11 * j11 + j12 * j12);
    std::println("cov2d  a={:.3f} b={:.3f} c={:.3f}", a, b, c);
    a += 0.3f;                                     // low-pass dilation
    c += 0.3f;
    const float det = a * c - b * b;
    std::println("conic  a={:.5f} b={:.5f} c={:.5f}", c / det, -b / det, a / det);
    const float mid = 0.5f * (a + c);
    const float l1 = mid + std::sqrt(std::max(0.1f, mid * mid - det));
    std::println("radius {} px", static_cast<int>(std::ceil(3.0f * std::sqrt(l1))));
}

The splat sits off-axis, so J's third column couples depth into x and y: the footprint is slightly elongated and tilted even though the 3D Gaussian is round.

4. Tiles: who touches what

The screen is cut into 16×16 tiles. Each Gaussian's 3σ square — centre m, half-size r — is converted to a rectangle of tiles (getRect):

min = ⌊(m − r) / 16⌋, max = ⌊(m + r + 15) / 16⌋, each clamped to [0, tiles],

and covers tiles [min, max). An empty rectangle (off-screen, or behind the near plane z ≤ 0.2) means the Gaussian is culled. tiles_touched = (max.x − min.x)·(max.y − min.y).

5. One sort for every tile

Now each (Gaussian, tile) pair becomes a 64-bit key: tile id in the high 32 bits, depth in the low 32. Sorting all keys at once groups them by tile and orders each tile's list by depth.

The reference then (1) prefix-sums tiles_touched so each Gaussian knows where to write its duplicates, (2) writes them, (3) runs CUB's DeviceRadixSort::SortPairs over only 32 + bit_width(tiles) key bits, and (4) finds each tile's range [start, end) by marking where the high word changes between neighbours.

#include <algorithm>
#include <bit>
#include <cstdint>
#include <print>
#include <utility>
#include <vector>

// Order-preserving for every finite float (the reference only needs depth > 0).
std::uint32_t depth_key(float d) {
    const std::uint32_t b = std::bit_cast<std::uint32_t>(d);
    return (b & 0x80000000u) ? ~b : (b | 0x80000000u);
}

int main() {
    struct Splat { float depth; int tile_lo, tile_hi; };   // tiles [lo, hi) of a 1-row grid
    const std::vector<Splat> splats{{3.0f, 0, 2}, {1.5f, 1, 2}, {4.2f, 0, 3}};
    std::vector<std::pair<std::uint64_t, std::uint32_t>> keys;
    for (std::uint32_t i = 0; i < splats.size(); ++i)
        for (int t = splats[i].tile_lo; t < splats[i].tile_hi; ++t)
            keys.emplace_back(std::uint64_t(t) << 32 | depth_key(splats[i].depth), i);
    std::ranges::sort(keys);
    for (std::size_t k = 0; k < keys.size(); ++k) {
        const auto tile = keys[k].first >> 32;
        const bool starts = k == 0 || (keys[k - 1].first >> 32) != tile;
        std::println("{} tile {} splat {}{}", k, tile, char('A' + keys[k].second), starts ? "  <- range start" : "");
    }
}

Note std::uint64_t(t) << 32: shifting a 32-bit value by 32 is undefined behaviour, a classic bug when porting key <<= 32 code that relied on key already being 64-bit.

6. Blending: front to back, per pixel

Each pixel walks its tile's range — nearest splat first — and composites with the front-to-back form of the volume rendering sum (gfx-volume-rendering):

C = Σᵢ cᵢ αᵢ Tᵢ + T_N · background, Tᵢ = Π_{j<i} (1 − αⱼ),

where αᵢ is the splat's opacity times its 2D Gaussian at the pixel. With d = mean2d − pixel and the conic (a, b, c), renderCUDA computes:

float power = -0.5f * (con.x * d.x * d.x + con.z * d.y * d.y) - con.y * d.x * d.y;
if (power > 0.0f) continue;                       // numerically impossible for a PD conic; guard
float alpha = min(0.99f, con_o.w * exp(power));   // cap: T never reaches 0 (the backward divides by 1 - α)
if (alpha < 1.0f / 255.0f) continue;              // invisible in 8-bit output: skip
float test_T = T * (1 - alpha);
if (test_T < 0.0001f) { done = true; continue; }  // opaque enough: stop, WITHOUT adding this splat
for (int ch = 0; ch < C; ch++) C[ch] += features[ch] * alpha * T;
T = test_T;

and finally out = C + T * bg. The order of those last lines is load-bearing: splat i is weighted by Tᵢ, the transmittance in front of it, so the update of T comes last. Pixel centres sit at integer coordinates: pixel (i, j) samples (i, j).

#include <algorithm>
#include <cmath>
#include <print>

struct Splat2D { float mx, my, ca, cb, cc, opacity, r, g, b; };   // conic (ca, cb, cc)

int main() {
    // Two splats over pixel (10, 10), nearest first, white background.
    const Splat2D splats[] = {{10.5f, 10.0f, 0.25f, 0.0f, 0.25f, 0.6f, 1, 0, 0},
                              {11.0f, 9.0f, 0.1f, 0.02f, 0.1f, 0.9f, 0, 0, 1}};
    const float px = 10, py = 10;
    float T = 1, C[3] = {0, 0, 0};
    for (const Splat2D& s : splats) {
        const float dx = s.mx - px, dy = s.my - py;
        const float power = -0.5f * (s.ca * dx * dx + s.cc * dy * dy) - s.cb * dx * dy;
        if (power > 0) continue;
        const float alpha = std::min(0.99f, s.opacity * std::exp(power));
        if (alpha < 1.0f / 255.0f) continue;
        const float test_T = T * (1 - alpha);
        if (test_T < 1e-4f) break;
        C[0] += s.r * alpha * T;
        C[1] += s.g * alpha * T;
        C[2] += s.b * alpha * T;
        std::println("alpha {:.4f}  weight {:.4f}", alpha, alpha * T);
        T = test_T;
    }
    std::println("pixel ({:.4f}, {:.4f}, {:.4f})  T = {:.4f}", C[0] + T, C[1] + T, C[2] + T, T);
}

On the GPU all 256 threads of a tile load the range in batches of 256 splats into shared memory together, and a block retires when every pixel is done. The reference also records, per pixel, the index of the last contributing splat and the final T: the backward pass needs both.

7. What's approximate, and what came after

• Popping. Sorting by the depth of each splat's centre, once per tile, is an approximation of the true order along each ray. When two overlapping splats' centre depths cross as the camera moves, their order flips for the whole tile at once and the colour jumps. StopThePop (Radl et al., SIGGRAPH 2024) sorts per pixel with a hierarchical, mostly-in-registers sort; 3D Gaussian Ray Tracing (3DGRT) and 3DGUT evaluate Gaussians along each ray in order.
• Loose bounds. The 3σ square wastes work on long, thin, diagonal splats. gsplat and follow-ups tighten it (per-axis extents, opacity-aware cut-offs: a splat with o = 0.05 is below 1/255 well before 3σ, tile–ellipse intersection tests).
• The linearised projection. J is exact only at the mean; wide-angle and fisheye cameras need better (gsplat's unscented-transform mode, 3DGUT).

gsplat vs the reference. Same math, different packaging: gsplat exposes fully_fused_projection + isect_tiles + rasterize_to_pixels, takes view matrices and intrinsics K instead of FoV + projection matrix, has an eps2d parameter for the 0.3, caps α at 0.999 instead of 0.99, offers a "packed" sparse mode, and computes exact gradients where the reference takes shortcuts (next skill).

In Python

import torch
from gsplat import rasterization

colors, alphas, meta = rasterization(
    means, quats, torch.exp(log_scales), torch.sigmoid(opacity_logits),
    sh_coeffs, viewmats, Ks, width, height, sh_degree=3,
    eps2d=0.3, rasterize_mode="classic")       # "antialiased" = Mip-Splatting's 2D filter

In an interview

"Walk me through rendering a Gaussian splat" — be ready to say: Σ = R S Sᵀ Rᵀ and why; Σ′ = J W Σ Wᵀ Jᵀ and what J is; the 0.3 dilation; conic and 3σ radius from eigenvalues; duplicate per tile, one radix sort on (tile, depth) keys; front-to-back with α = min(0.99, o·G), early termination at T < 10⁻⁴; and that per-tile sorting causes popping. Then the follow-up: what would you change for a mobile viewer? (fewer SH bands, half floats, sort on a worker, bounded splat count, compression).