Compressing and streaming Gaussian splats

A trained 3D Gaussian Splatting scene (Kerbl et al. 2023) is a list of millions of anisotropic Gaussians, each with a position, a shape, an opacity and a view-dependent colour. Rendering it is cheap; shipping it is not. The reference .ply of a single room is half a gigabyte, and a company whose product is "photorealistic 3D on every device" has to get that onto a phone over a mobile connection, fast, and start drawing before it has all arrived.

This lesson builds the toolbox: where the bytes go, IEEE half floats done right, quantisation against local bounds (with Morton order making the bounds local), the "smallest three" trick for rotations, log and logit spaces, spherical-harmonic pruning, vector quantisation with k-means codebooks, entropy coding, how to measure what you lost — and level-of-detail streaming.

1. Where the bytes go

The reference implementation writes one .ply vertex per splat with 62 float32 properties:

Drop the normals and you have 59 floats = 236 bytes. Three million splats — an ordinary captured room — is 708 MB. The obvious observations drive every compressed format:

• SH is the elephant: 45 of 59 floats. Dropping degrees 1–3 entirely (diffuse colour only) already cuts 76 %; most compressed formats keep them as an optional, heavily quantised stream.
• float32 is wildly over-precise. Positions need ~1 mm on a 10 m scene, colours 8 bits, rotations a fraction of a degree.
• Neighbouring splats are similar. Sort them well and generic compressors find long runs of similar bytes.

#include <print>

int main() {
    const double splats = 3e6;
    const double raw = 59 * 4;             // floats → bytes, normals dropped
    const double no_sh = (59 - 45) * 4;    // keep only SH degree 0
    const double packed = 16 + 48.0 / 256; // the chunked format of this lesson
    std::println("raw     {:6.1f} MB", splats * raw / 1e6);
    std::println("no SH   {:6.1f} MB", splats * no_sh / 1e6);
    std::println("packed  {:6.1f} MB  ({:.1f}x smaller than raw)", splats * packed / 1e6, raw / packed);
}

2. Half floats, rounded correctly

The first, easy win is fp16 (IEEE binary16): 1 sign bit, 5 exponent bits (bias 15), 10 mantissa bits. It holds about 3.3 significant decimal digits, normal numbers from 2⁻¹⁴ ≈ 6.1·10⁻⁵ to 65504, subnormals down to 2⁻²⁴ ≈ 6·10⁻⁸. GPUs read it natively, which makes it the standard format for SH coefficients and other "precise enough" attributes.

The conversion has four regimes, keyed by the re-biased exponent e = exp − 127 + 15:

Rounding is the part people get wrong. IEEE's default is round to nearest, ties to even: look at the dropped bits; above half an ulp round up, below round down, exactly half round so the last kept bit is 0. Plain truncation (mant >> 13) is wrong for about half of all inputs and biases every value toward zero — millions of slightly-too-dark colours. Two details make the code short: the rounding increment is applied to the packed exponent|mantissa value, so a mantissa overflow carries into the exponent (1.111…1 rounds up to 2.0), and 0x7BFF + 1 = 0x7C00 is exactly infinity, so 65520 overflows correctly.

#include <bit>
#include <cassert>
#include <cstdint>
#include <print>

std::uint16_t float_to_half(float f) {
    const std::uint32_t x = std::bit_cast<std::uint32_t>(f);
    const std::uint32_t sign = (x >> 16) & 0x8000u;
    const std::uint32_t exp = (x >> 23) & 0xFFu;
    std::uint32_t mant = x & 0x7FFFFFu;
    if (exp == 0xFF) return sign | 0x7C00u | (mant ? 0x200u | (mant >> 13) : 0u);
    const int e = int(exp) - 127 + 15;
    if (e >= 31) return sign | 0x7C00u;
    if (e <= 0) {                                    // subnormal half
        if (e < -10) return sign;
        mant |= 0x800000u;
        const int shift = 14 - e;
        std::uint32_t h = mant >> shift;
        const std::uint32_t rem = mant & ((1u << shift) - 1), half = 1u << (shift - 1);
        if (rem > half || (rem == half && (h & 1))) ++h;
        return sign | h;
    }
    std::uint32_t h = (std::uint32_t(e) << 10) | (mant >> 13);
    const std::uint32_t rem = mant & 0x1FFFu;
    if (rem > 0x1000u || (rem == 0x1000u && (h & 1))) ++h;   // ties to even
    return sign | h;
}

int main() {
    assert(float_to_half(1.0f) == 0x3C00);
    assert(float_to_half(0.1f) == 0x2E66);
    assert(float_to_half(2049.0f) == 0x6800);        // tie → even (2048)
    assert(float_to_half(2051.0f) == 0x6802);        // tie → even (2052)
    assert(float_to_half(65520.0f) == 0x7C00);       // rounds up to infinity
    assert(float_to_half(0x1p-25f) == 0x0000);       // tie between 0 and 2^-24 → 0
    assert(float_to_half(0x1.8p-25f) == 0x0001);
    std::println("all fp16 checks pass");
}

In C++23 you can also write static_cast<std::float16_t>(f) (<stdfloat>) or the _Float16 extension and let the compiler (or the F16C instruction vcvtps2ph) do it — that is exactly what this skill's tests use as an oracle. Writing it yourself is still a common interview question, and you need it on targets without the type.

3. Quantise against local bounds

fp16 has a fatal flaw for positions: its precision is relative. At x = 1000 the spacing between halves is 0.5 — useless. What you want is a fixed number of steps across a box:

q = round(t · (2^bits − 1)),  t = clamp((v − lo) / (hi − lo), 0, 1)
v' = lo + (hi − lo) · q / (2^bits − 1)        |v − v'| ≤ (hi − lo) / (2 · (2^bits − 1))

Three things can go wrong in those two lines, and all three are classic:

• hi == lo (every splat in the box has the same y) makes t = 0/0 = NaN, and converting NaN to an integer is undefined behaviour. Check !(hi > lo) first.
• Without the clamp, v slightly outside the range (colours from SH often are) produces 2^bits or more — one bit too many, which silently corrupts the neighbouring field of a packed word.
• 1u << bits is UB for bits == 32.

The error bound is set by the box. One box for a whole 100 m city scene gives 11-bit steps of 5 cm — visibly jittery. So we cut the scene into small boxes. Which splats go into a box? If you take them in file order, 256 random splats span the whole scene and nothing is gained. Sort them along a Morton (Z-order) curve first (see gfx-morton-octree): interleave the bits of the quantised x, y, z, sort by that code, and consecutive splats become spatial neighbours. Then every run of 256 has a tight bounding box.

This is the PlayCanvas "compressed .ply" used by SuperSplat: chunks of 256 Morton-sorted splats; per chunk a 48-byte header (min/max of position and of log-scale, 12 floats); per splat 16 bytes — position as one uint32 split 11-10-11 bits, rotation 2+10+10+10, log-scale 11-10-11, colour + opacity as RGBA8. The header costs 0.19 bytes per splat; the result is about 15× smaller than the full .ply (≈4× smaller than a .ply with SH stripped).

Why 11-10-11 rather than 10-10-10? It uses all 32 bits; giving the odd bit to x and z (the horizontal axes of most captured scenes, which are wider than they are tall) is a heuristic. Watch the effect of local bounds:

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

float roundtrip(float v, float lo, float hi, int bits) {
    const float maxq = float((1u << bits) - 1u);
    if (!(hi > lo)) return lo;
    const float t = std::clamp((v - lo) / (hi - lo), 0.0f, 1.0f);
    const auto q = std::uint32_t(t * maxq + 0.5f);
    return lo + (hi - lo) * (float(q) / maxq);
}

int main() {
    const float x = 37.4131f;   // a splat 37 m from the origin
    std::println("global box [-50, 50], 11 bits: error {:.2e} m", std::fabs(roundtrip(x, -50, 50, 11) - x));
    std::println("chunk  box [37, 38],  11 bits: error {:.2e} m", std::fabs(roundtrip(x, 37, 38, 11) - x));
    std::println("flat chunk [2, 2]:            error {:.2e} m", std::fabs(roundtrip(2.0f, 2, 2, 11) - 2.0f));
}

Note the bounds themselves stay float32: rounding a bound to fp16 could move it inward, and then the extreme splats no longer fit.

4. Rotations in 32 bits: "smallest three"

A unit quaternion (w, x, y, z) has three degrees of freedom; the fourth component is ±√(1 − others²). Which one to drop? The largest in magnitude, for two reasons: its value is the most precisely recoverable (the square root is flattest near 1), and the other three are then bounded: if |c| > 1/√2 for some non-largest c, the largest would also exceed 1/√2 and the squared sum would exceed 1. So each stored component lies in [−1/√2, 1/√2], and 10 bits over that range give a step of 0.0014. Glenn Fiedler popularised this for networked physics ("Snapshot Compression"); SPZ v3+ and the PlayCanvas format use it.

The decoder's square root is always positive. If the dropped component was negative, the encoder must negate the whole quaternion first — q and −q are the same rotation (the double cover of SO(3)). Negating only the dropped component gives a different rotation; forgetting it entirely flips the sign of one component of four. Both pass a test that only uses quaternions near identity.

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

using Quat = std::array<float, 4>;
constexpr float k = std::numbers::sqrt2_v<float> / 2;

std::uint32_t q10(float v) { return std::uint32_t((std::clamp(v, -k, k) + k) / (2 * k) * 1023.0f + 0.5f); }
float d10(std::uint32_t q) { return -k + 2 * k * (float(q) / 1023.0f); }

std::uint32_t pack(Quat q) {
    int L = 0;
    for (int i = 1; i < 4; ++i) if (std::fabs(q[i]) > std::fabs(q[L])) L = i;
    if (q[L] < 0) for (float& c : q) c = -c;          // −q is the same rotation
    std::uint32_t p = std::uint32_t(L) << 30;
    for (int i = 0, s = 20; i < 4; ++i) if (i != L) { p |= q10(q[i]) << s; s -= 10; }
    return p;
}

Quat unpack(std::uint32_t p) {
    const int L = int(p >> 30);
    Quat q{};
    float sum = 0;
    for (int i = 0, s = 20; i < 4; ++i) if (i != L) { q[i] = d10((p >> s) & 1023u); sum += q[i] * q[i]; s -= 10; }
    q[L] = std::sqrt(std::max(0.0f, 1 - sum));
    return q;
}

int main() {
    std::uint32_t seed = 1;
    auto rnd = [&] { seed = seed * 1664525u + 1013904223u; return float(seed >> 8) / 8388608.0f - 1.0f; };
    double worst = 0;
    for (int n = 0; n < 100000; ++n) {
        Quat q{rnd(), rnd(), rnd(), rnd()};
        const float len = std::sqrt(q[0] * q[0] + q[1] * q[1] + q[2] * q[2] + q[3] * q[3]);
        if (len < 1e-3f) continue;
        for (float& c : q) c /= len;
        const Quat r = unpack(pack(q));
        double dot = 0;
        for (int i = 0; i < 4; ++i) dot += double(q[i]) * r[i];
        worst = std::max(worst, 2 * std::acos(std::min(1.0, std::fabs(dot))) * 180 / std::numbers::pi);
    }
    std::println("worst rotation error over 100k quaternions: {:.3f} degrees", worst);
}

5. Opacity, scale and colour

The 3DGS optimiser stores each attribute in the space where gradient descent behaves, and those spaces happen to be good for quantisation too:

• Scale is stored as log(σ): splats range from sub-millimetre to metres, and a uniform step in log space is a constant relative error, e^{Δ/2} − 1. Quantise log σ per chunk (11-10-11 bits) and exp on decode.
• Opacity is stored as a logit a = ln(α / (1 − α)) and passed through a sigmoid. For storage, apply the sigmoid and keep α in 8 bits — the eye cares about α linearly — or quantise the logit when you need precision near 0 and 1 (nearly transparent splats matter: there are many of them).
• Colour: degree-0 SH f_dc maps to RGB as 0.5 + C₀ · f_dc with C₀ = 1/(2√π) ≈ 0.2821. Store 8 bits per channel after clamping — values outside [0, 1] occur.
• Higher-order SH (degrees 1–3, see gfx-spherical-harmonics): the coefficients are small and their errors are seen only as subtle view-dependent shading. Options, roughly in order of aggressiveness: fp16; 8 bits on [−1, 1]; SPZ's 5 bits for degree 1 and 4 bits for degrees 2–3; distil to a lower degree (LightGaussian drops degree 3, re-training so degree 2 absorbs it); drop entirely for splats whose SH energy is tiny; or a codebook (next section).

6. Vector quantisation: k-means codebooks

Scalar quantisation treats each number alone. Many splats, though, have the same kind of view-dependent colour. Vector quantisation clusters the 45-dimensional SH vectors with k-means into K centroids (K = 4096, say), stores the codebook once, and gives each splat just a 12-bit index: 2 million × 45 × 4 B = 360 MB becomes 3 MB of indices plus a 369 KB codebook. CompGS (first titled Compact3D, Navaneet et al. 2023) and LightGaussian (Fan et al. 2023) do this for colour and for the covariance (scale + rotation), then fine-tune the codebook through the renderer so the errors that matter shrink; Lee et al. (CVPR 2024) use residual VQ for geometry and a hash grid for colour.

k-means itself is Lloyd's algorithm: assign each vector to its nearest centroid, move each centroid to the mean of its members, repeat. The cost is the assignment, O(n · K · d) per iteration — at 2M × 4096 × 45 that is 3.7·10¹¹ multiply-adds, so production code runs it on the GPU, on a random subsample, or with a kd-tree/product-quantised search for the nearest centroid.

7. Entropy coding and layout

Quantisation makes values small; entropy coding makes the stream short. Generic compressors (gzip/DEFLATE, zstd) find repeated byte patterns, so help them:

• Columnar layout. Write all positions, then all rotations, then all colours (structure of arrays), not 16-byte records. Each column is homogeneous, and SPZ v4 compresses each column as its own zstd stream (in parallel).
• Spatial order. After a Morton sort, neighbouring splats have similar colours, scales and quantised positions; delta-encoding positions within a chunk turns them into small numbers.
• Learned / context models. Research codecs (HAC, Chen et al. 2024) predict each attribute's distribution from a hash-grid context and arithmetic-code it; Self-Organizing Gaussians (Morgenstern et al.) sort splats into 2D grids so standard image codecs compress them.

SPZ (Niantic, open source) is the real-world reference point: 24-bit fixed-point positions (fractional bits in the header), 8-bit log scales, smallest-three rotations with 10 bits per component (v3+; v2 stored x, y, z in 8 bits and derived w), 8-bit alpha and colour, SH as 8-bit values keeping 5 bits for degree 1 and 4 for higher degrees, columns compressed separately — about 10× smaller than .ply "with virtually no perceptible loss".

8. Measuring what you lost

Compression papers report three image metrics on held-out views:

• MSE: mean of squared pixel differences.
• PSNR = 10 · log10(MAX² / MSE) dB. MAX is the largest value the data can take: 255 for 8-bit images, 1 for images normalised to [0, 1]. Mix them up and every number moves by 20 · log10 255 ≈ 48.13 dB. Each halving of RMS error adds 6.02 dB; good 3DGS reconstructions sit around 25–35 dB on standard benchmarks, and a codec is judged by how many dB it costs (typically < 0.5).
• SSIM (Wang et al. 2004) compares local means, variances and covariance in small windows — structure rather than raw error — and is 1 for identical images. The 3DGS training loss is 0.8·L1 + 0.2·(1 − SSIM). LPIPS (Zhang et al. 2018) compares deep features and tracks perceived quality best.

#include <cmath>
#include <cstdint>
#include <limits>
#include <print>
#include <span>
#include <vector>

double psnr_u8(std::span<const std::uint8_t> a, std::span<const std::uint8_t> b) {
    double sum = 0;
    for (std::size_t i = 0; i < a.size(); ++i) {
        const double d = double(a[i]) - double(b[i]);
        sum += d * d;
    }
    const double mse = sum / double(a.size());
    if (mse == 0) return std::numeric_limits<double>::infinity();
    return 10.0 * std::log10(255.0 * 255.0 / mse);
}

int main() {
    std::vector<std::uint8_t> ref(256), q4(256);
    for (int i = 0; i < 256; ++i) {
        ref[i] = std::uint8_t(i);
        q4[i] = std::uint8_t((i >> 4) * 17);           // keep 4 bits, expand back to 0..255
    }
    std::println("8 -> 4 bits per channel: PSNR {:.2f} dB", psnr_u8(ref, q4));
    std::println("identical:               PSNR {}", psnr_u8(ref, ref));
}

9. Level of detail and progressive streaming

The final trick isn't about bytes per splat but about which splats first. A phone rendering at 1080p cannot resolve millions of sub-pixel splats in the distance anyway. Build a hierarchy: an octree over the scene (see gfx-morton-octree) whose inner nodes hold merged splats — one Gaussian fitted to its children's mixture (moment matching: weighted mean of positions, covariance of the mixture). Then:

• send the top levels first: a few thousand coarse splats cover the whole view in the first hundred kilobytes;
• refine nodes by screen-space error — expand a node when its splats would cover more than a pixel or two — so nearby detail arrives before distant detail;
• render whatever has arrived; each chunk replaces its parent's coarse splats.

Hierarchical 3DGS (Kerbl et al. 2024) and Octree-GS (Ren et al. 2024) build such hierarchies for city-scale scenes; streaming viewers do the same with Morton-sorted chunks as the unit of transfer.

In an interview

• "How would you compress this .ply?" Start from the byte budget (SH is 76 %), then: Morton sort → chunks → per-chunk quantisation; smallest-three rotations; log scales, sigmoid opacity, 8-bit colour; SH degree reduction or codebooks; columnar layout + zstd. Mention PSNR/SSIM/LPIPS before and after.
• "Write float → half." Know the four regimes and round-to-nearest-even; test 65520 → ∞, 2049 → 2048, 2⁻²⁵ → 0.
• "Why drop the largest quaternion component?" Bounded range for the rest (±1/√2), most stable reconstruction, and remember the sign flip.
• Numerical traps to name unprompted: zero-extent ranges, clamping before packing, NaN-to-int conversion, 1 << 32, PSNR's MAX.

In Python (NumPy)

import numpy as np
h = np.float32([1.0, 2049.0, 65520.0]).astype(np.float16)   # round to nearest even
print(h.view(np.uint16))                                     # [15360 26624 31744]

lo, hi = pos.min(0), pos.max(0)                               # per chunk, shape (3,)
ext = np.where(hi > lo, hi - lo, 1.0)                         # guard zero extent
q = np.rint(np.clip((pos - lo) / ext, 0, 1) * [2047, 1023, 2047]).astype(np.uint32)
packed = q[:, 0] << 21 | q[:, 1] << 11 | q[:, 2]

PyTorch training code (gsplat, the reference diff-gaussian-rasterization) keeps the raw parameters in float32 and applies exp, sigmoid and quaternion normalisation in the forward pass — which is why the .ply stores log scales, logit opacities and unnormalised quaternions.