Gaussian splat files with NumPy

A trained 3D Gaussian splatting scene is nothing but a long list of Gaussians, each with a position, a shape, an opacity and a colour that changes with the viewing direction. The way that list reaches a viewer — SuperSplat, the SIBR viewer, a web player, a streaming server — is almost always the .ply file written by GaussianModel.save_ply in graphdeco-inria/gaussian-splatting. There is no other spec: that function is the format. This lesson reads and writes it with NumPy — one structured dtype and np.frombuffer, no loops — and covers every quirk that trips loaders: the header's byte offset, typed properties, the channel-major SH coefficients, the activations that must be applied on load, and the quaternion order.

1. What is stored per splat

For each Gaussian the optimiser holds these tensors (scene/gaussian_model.py):

plus nx ny nz, three normals that are always written as zero (the code started from a point-cloud writer). The file stores the raw optimiser parameters, before activation: Adam updates unconstrained reals, and sigmoid/exp map them to valid opacities and scales. That is the single most important fact about the format, and section 5 is about it.

2. The PLY header

PLY (Greg Turk's Stanford format) is an ASCII header followed by a body. The header declares elements (here just vertex, sometimes followed by an empty face element) and their properties, each with a type:

ply
format binary_little_endian 1.0
comment written by 3DGS save_ply (via plyfile)
element vertex 1843209
property float x
property float y
property float z
property float nx
...
property float rot_3
end_header
<binary records>

The body of a binary_little_endian file is count records, each the element's properties packed back to back in declaration order, without padding. So the header alone tells you the byte layout. Parsing it is a loop over lines:

PLY_TYPES = {"char": "i1", "uchar": "u1", "short": "i2", "ushort": "u2", "int": "i4",
             "uint": "u4", "float": "f4", "double": "f8", "float32": "f4", "uint8": "u1"}

def parse_header(data):
    end = data.find(b"end_header\n")
    if not data.startswith(b"ply\n") or end < 0:
        raise ValueError("not a PLY file")
    count, props, element = None, [], None
    for line in data[:end].decode("ascii").split("\n")[1:]:
        w = line.split()
        if not w or w[0] in ("comment", "obj_info"):
            continue
        if w[0] == "format" and w[1] != "binary_little_endian":
            raise ValueError(w[1])
        if w[0] == "element":
            element = w[1]
            if element == "vertex":
                count = int(w[2])
        elif w[0] == "property" and element == "vertex":
            props.append((w[2], "<" + PLY_TYPES[w[1]]))
    return count, props, end + len(b"end_header\n")

header = (b"ply\nformat binary_little_endian 1.0\ncomment hi\nelement vertex 2\n"
          b"property float x\nproperty uchar red\nproperty double opacity\n"
          b"element face 0\nproperty list uchar int vertex_indices\nend_header\n")
print(parse_header(header))
print("header is", len(header), "bytes")

Three things to get exactly right. The body starts after the newline that ends end_header — stop one byte short and every float is assembled from bytes of two neighbours. Property types are declared, not assumed: files from 3DGS are all float, but exports from other tools add uchar red green blue or double coordinates, and a loader that assumes 4 bytes per property reads everything after the first such column from the wrong offset. And the byte order comes from the format line, not from your machine.

3. Structured dtypes and np.frombuffer

A NumPy structured dtype describes one record: a list of (name, type) fields, packed like a C struct with #pragma pack(1). Its itemsize is the record size, and np.frombuffer turns the body into an array of records without copying or looping:

import numpy as np

props = [("x", "<f4"), ("red", "u1"), ("opacity", "<f8")]
dt = np.dtype(props)
print("itemsize", dt.itemsize, "offsets", [dt.fields[n][1] for n in dt.names])

body = np.array([(1.5, 200, -2.0), (-3.0, 7, 4.5)], dtype=dt).tobytes()
v = np.frombuffer(body, dtype=dt, count=2, offset=0)
print(v["x"], v["red"], v["opacity"])
print("1.0f little-endian bytes:", np.float32(1.0).astype("<f4").tobytes().hex(),
      "read as big-endian:", np.frombuffer(np.float32(1.0).astype("<f4").tobytes(), ">f4")[0])

Field access v["x"] is a strided view: no copy, stride = itemsize. Two consequences. Views of a bytes object are read-only (bytes are immutable), so a loader must copy before anyone modifies a parameter — astype(np.float32) copies by default, but np.asarray(v["x"], np.float32) does not when the dtype already matches. And a field that sits at an odd offset (after a uchar) is unaligned; NumPy handles it, but copying into contiguous float32 arrays before heavy math is faster anyway.

4. Spherical harmonics and f_rest's channel-major order

Each splat's colour is a function of the viewing direction, stored as spherical-harmonic coefficients: K = (ℓ + 1)² per colour channel for degree ℓ — 1, 4, 9 or 16. Coefficient 0 (the "DC" term) is the view-independent base colour; the other K − 1 add view dependence. In memory 3DGS keeps them as (N, K, 3): coefficient-major, channel last. But save_ply writes

f_rest = self._features_rest.detach().transpose(1, 2).flatten(start_dim=1).cpu().numpy()

— transposed to (N, 3, K − 1) first. On disk the higher bands are channel-major: f_rest_0 … f_rest_{K−2} are all of red's coefficients, then green's, then blue's. load_ply undoes it with reshape((N, 3, K − 1)) followed by transpose(1, 2).

import numpy as np

N, K = 2, 4                                         # degree 1
rng = np.random.default_rng(0)
sh = rng.normal(size=(N, K, 3)).astype(np.float32)  # in-memory layout (N, K, 3)

f_rest = sh[:, 1:, :].transpose(0, 2, 1).reshape(N, -1)       # what save_ply writes
back = f_rest.reshape(N, 3, K - 1).transpose(0, 2, 1)         # what load_ply does
print("file order == channel-major:", np.array_equal(f_rest[0, :3], sh[0, 1:, 0]))
print("round trip exact:", np.array_equal(back, sh[:, 1:, :]))
print("wrong reshape equal?", np.array_equal(f_rest.reshape(N, K - 1, 3), sh[:, 1:, :]))

names = [f"f_rest_{i}" for i in range(45)]
print("sorted() order:", sorted(names)[:4], "...")
print("numeric order: ", sorted(names, key=lambda n: int(n.split("_")[-1]))[:4], "...")

That last pair is the second classic f_rest bug. 3DGS's load_ply sorts the f_rest property names — robust to exporters that reorder them — with key=lambda x: int(x.split('_')[-1]). Plain sorted() is lexicographic and puts f_rest_10 before f_rest_2; degree-1 files (f_rest_0..8) never show it. The number of f_rest fields also tells you the degree: 0, 9, 24, 45 for ℓ = 0…3; anything else is a broken file.

5. Activations on load

What the file holds is what the optimiser held; what the renderer needs is the activated value:

• opacity = sigmoid(logit). 3DGS initialises it with inverse_sigmoid(0.1) ≈ −2.2; trained values span roughly −10 to +10, and pruning pushes some far negative. The naive 1/(1 + exp(−x)) evaluates exp(+large) for x ≪ 0 — an overflow (a warning, or an error under np.errstate(over="raise")) even though the final 0 is right. The stable form exponentiates only −|x|.
• scales = exp(log-scale): standard deviations along the Gaussian's axes, in world units. Initialised to log(mean distance to the 3 nearest points).
• rotation: normalise the quaternion, which is w-first: 3DGS's build_rotation reads r = q[:, 0] as the real part, and gsplat uses wxyz as well. scipy's Rotation.from_quat wants xyzw — convert at that boundary only. A zero quaternion (a corrupted splat) has no direction; map it to the identity.
• colour: the view-independent part is 0.5 + C0·f_dc with C0 = 1/(2√π) = 0.28209479…, the value of the constant SH basis function. The 0.5 makes zero coefficients grey (3DGS's RGB2SH(rgb) = (rgb − 0.5)/C0 initialises f_dc from the point cloud's colours). The rasterizer clamps the result at 0 but not at 1:

// diff-gaussian-rasterization, computeColorFromSH (forward.cu)
glm::vec3 result = SH_C0 * sh[0];
if (deg > 0) { /* + bands 1..3, evaluated at the view direction */ }
result += 0.5f;
return glm::max(result, 0.0f);

import numpy as np

SH_C0 = 0.28209479177387814

def sigmoid(x):
    e = np.exp(-np.abs(x))                    # never exp of a positive number
    return np.where(x >= 0, 1 / (1 + e), e / (1 + e)).astype(x.dtype)

logits = np.array([-1000, -2.197, 0, 3, 1000], np.float32)
with np.errstate(over="raise"):
    print("opacity", [round(float(v), 4) for v in sigmoid(logits)])
    try:
        1 / (1 + np.exp(-logits))
    except FloatingPointError as e:
        print("naive sigmoid:", e)

print("scales ", [round(float(v), 4) for v in np.exp(np.array([0.0, np.log(2), -4.6], np.float32))])
q = np.array([[2, 0, 0, 0], [1, 1, 1, 1], [0, 0, 0, 0]], np.float32)
n = np.linalg.norm(q, axis=-1, keepdims=True)
print("quats  ", np.where(n > 0, q / np.where(n > 0, n, 1), [1, 0, 0, 0]).tolist())
print("rgb    ", np.maximum(0.5 + SH_C0 * np.array([[0.0, 1.0, -3.0]]), 0).round(4).tolist())

Note np.where evaluates both branches before selecting — so both must be safe on every element. That's why the sigmoid computes one e = exp(−|x|) and the quaternion code divides by where(n > 0, n, 1) rather than by n.

6. Writing and round-tripping

Writing is the same steps backwards: put the columns in the 3DGS order, write the header with property float <name> lines, and append the records. A (N, 62) float32 array in C order is the body — row-major means record after record:

import numpy as np

def write_ply(means, sh, opacity, scales, quats):
    n, k = sh.shape[:2]
    rest = sh[:, 1:, :].transpose(0, 2, 1).reshape(n, 3 * (k - 1))   # channel-major
    cols = np.concatenate([means, np.zeros((n, 3)), sh[:, 0, :], rest,
                           opacity[:, None], scales, quats], axis=1).astype("<f4")
    names = ["x", "y", "z", "nx", "ny", "nz"] + [f"f_dc_{i}" for i in range(3)]
    names += [f"f_rest_{i}" for i in range(rest.shape[1])] + ["opacity"]
    names += [f"scale_{i}" for i in range(3)] + [f"rot_{i}" for i in range(4)]
    head = f"ply\nformat binary_little_endian 1.0\nelement vertex {n}\n"
    head += "".join(f"property float {m}\n" for m in names) + "end_header\n"
    return head.encode("ascii") + cols.tobytes()

rng = np.random.default_rng(1)
n = 1000
data = write_ply(rng.normal(size=(n, 3)), rng.normal(size=(n, 16, 3)), rng.normal(size=n),
                 rng.normal(size=(n, 3)), rng.normal(size=(n, 4)))
body = data.index(b"end_header\n") + 11
print("header", body, "bytes; body", len(data) - body, "bytes =", (len(data) - body) // n, "per splat")

Test a writer against the layout, not only against your own reader: a reader and writer that both skip the f_rest transpose round-trip perfectly and still produce files every other viewer renders wrongly.

7. Sizes, and why everyone compresses

Bytes per splat at each SH degree, with normals:

A typical outdoor scene has 2–6 million splats: 0.5–1.5 GB of float32 at degree 3, far too much to stream to a phone. Compression schemes attack the table in order of size: drop the zero normals (−12 B), prune or codebook the higher SH bands (−180 B), then quantise the rest per spatial chunk — positions relative to a chunk's bounding box, log-scales and smallest-three quaternions in a few bits each (PlayCanvas's compressed .ply, Niantic's SPZ, the C++ compressing and streaming Gaussian splats skill). Even the antimatter15 .splat format — float32 position and scale, uint8 RGBA and quaternion, no SH: 32 bytes — is 7.75× smaller. Note that .splat stores activated, linear opacity and colour as bytes; feeding .ply logits to a .splat code path (or vice versa) is a common conversion bug.

8. Conventions to check before trusting a loader

The build task is the loader and writer from this lesson, made robust: typed properties in any order, extra fields and elements, degrees 0–3, overflow-free activations, writable copies, exact round trips — and the bytes-per-splat arithmetic that motivates compressing it all.