tsplat

Run Gaussian Splatting in your terminal, CPU only, even works over SSH

tsplat renders 3D Gaussian Splatting scenes directly in your terminal using Unicode half-blocks or any supported graphics protocol. written in rust and is CPU only for now, doesn't require GPU or display servers and it even works over SSH

Installation

Requires Rust

cargo install --git https://github.com/darshanmakwana412/tsplat

Or clone and build locally:

git clone https://github.com/darshanmakwana412/tsplat
cd tsplat
cargo build --release

Quick start

You need an INRIA 3DGS .ply scene. The pretrained garden scene is a good first test:

tsplat path/to/scene.ply

# Raise or remove the default 200 k splat cap
tsplat path/to/scene.ply --max-splats 500000
tsplat path/to/scene.ply --no-cap

# Scene looks washed-out when
tsplat path/to/scene.ply --raw-opacity

Note: If you terminal does not support any graphics protocol it will fallback to half block rendering which might look minecrafty

Controls

Key	Action
W A S D	Pan (camera-relative)
J / K or left/right arrows	Yaw left / right
H / L or up/down arrows	Pitch up / down
+ / - or mouse scroll	Zoom in / out
O	Toggle auto-orbit (smooth yaw loop, ideal for recordings)
Tab	Toggle HUD — arrows navigate rows and adjust values while open
q / Esc / Ctrl-C	Quit

Benchmarks

cargo run --release --bin bench_forward
cargo run --release --bin bench_forward -- --threads 1,2,4,8 --splats 200000
cargo run --release --bin bench_forward -- --ply path/to/scene.ply

---

Tutorial on Gaussian Splatting

I recently discovered some of my notes related to it and decided to digitize it this weekend, along the way I reimplemented the forward rasterization pass in rust and decided it would be fun to write a tutorial explaining gaussian splatting to everyone, so here it is

• What is a gaussian splat?
• The forward pass pipeline
• Step 1: Projecting Splats
• 1.1: building the 3D covariance matrix
• 1.2: transforming into view space
• 1.3: projecting to 2D
• Step 2: computing bounding boxes
• Step 3: depth sort
• Step 4: tile binning
• Step 5: alpha compositing
• tiled parallel compositing
• Optimizations
• fast approximate exp
• row-level early-out
• scratch buffer reuse

what is a gaussian splat?

a 3D Gaussian splat is an oriented ellipsoid in space that carries some color and opacity. you can think of it as a fuzzy colored blob. a scene is made of hundreds of thousands of these blobs, and when you look at them from a particular viewpoint, they overlap and blend to form the final image

We represent each gaussian with these attributes:

pub struct Splat {
    pub pos: Vec3,      // center position in world space
    pub scale: Vec3,    // size along each local axis
    pub rot: Quat,      // orientation as a unit quaternion
    pub color: Vec3,    // RGB color (already decoded from spherical harmonics)
    pub opacity: f32,   // how opaque this blob is, in [0, 1]
}

scales are stored in log space, opacities as logits, colors as spherical harmonics coefficients and quaternions are normzlied to unit length to ensure the values lie within their respective range

spherical harmonic (SH) coefficients are just a frequency-domain representation of a color function defined over the unit sphere, now why spherical harmonics? because in the real world, the color of a surface depends on the viewing direction. SH coefficients encode this view-dependent appearance compactly.

SH functions are organized in bands (like octaves in music), as you go higher up in the bands you have more coefficients and thus they capture more finer details, the INRIA 3DGS format stores up to band 3 (48 coefficients per splat for RGB)

To decode bash 0, the band-0 SH basis function is $Y_0^0 = \frac{1}{2\sqrt{\pi}} \approx 0.282$. the conversion from SH coefficient to RGB is:

$$\text{color} = \text{clamp}\left(0.5 + C_0 \cdot f_{dc},\ 0,\ 1\right)$$

where $C_0 = Y_0^0$ and $f_{dc}$ is the 3-component DC coefficient from the file.

pub const SH_C0: f32 = 0.28209479177387814;

pub fn sh_band0_to_rgb(f_dc: Vec3) -> Vec3 {
    (Vec3::splat(0.5) + SH_C0 * f_dc).clamp(Vec3::ZERO, Vec3::ONE)
}

the forward pass pipeline

the forward pass turns a list of 3D Gaussians + a camera into a 2D image. here is an overview of the rendering pipeline:

Step 1: Projecting Splats

1.1: building the 3D covariance matrix

for each splat given the raw (scale, rotation) pairs we need to construct a 3D covariance matrix $\Sigma$ that describes the shape and orientation of the Gaussian in world space. the formula is:

$$\Sigma = R \cdot S \cdot S^T \cdot R^T$$

where R is the 3×3 rotation matrix from the quaternion, and S is a diagonal matrix of scales. if we let M = R·S, this simplifies to:

$$\Sigma = M \cdot M^T$$

let r_mat = Mat3::from_quat(s.rot);
let s_mat = Mat3::from_diagonal(s.scale);
let m = r_mat * s_mat;
let cov3d = m * m.transpose();

---

Note: why dowe decompose the covariance this way?

Covariance matrices have physical meaning only when they are positive semi-definite. gradient descent cannot easily be constrained to produce valid matrices, by expressing the covariance as $M \cdot M^T$, it is guaranteed to be positive semi-definite, a matrix of the form $A^T A$ always is. this is a reparametrization trick: we optimize scale and rotation separately, which are unconstrained, and the covariance we derive from them is always valid

---

what does this matrix actually look like? for a splat with scale = (0.1, 0.05, 0.02) and identity rotation:

$$ \Sigma = \begin{pmatrix} 0.01 & 0 & 0 \ 0 & 0.0025 & 0 \ 0 & 0 & 0.0004 \end{pmatrix} \quad \begin{aligned} &= \text{diag}(0.1^2,\; 0.05^2,\; 0.02^2) \end{aligned} $$

with identity rotation, it is just the squared scales on the diagonal, an axis-aligned ellipsoid

1.2: transforming into view space

the 3D covariance we just computed lives in world space. to project it onto the camera's image plane, we first need to rotate it into view space, the coordinate system where the camera is at the origin, looking down −z

for the splat center, this is just a matrix-vector multiply with the 4×4 view matrix:

let p_view4 = view * Vec4::new(s.pos.x, s.pos.y, s.pos.z, 1.0);
let p_view = Vec3::new(p_view4.x, p_view4.y, p_view4.z);
if p_view.z > -znear || p_view.z < -zfar {
    return None;
}
let zc = -p_view.z;

note zc = -p_view.z. our view space is right-handed with the camera looking down −z, so points in front of the camera have negative z. we use zc (positive in front) as the depth for sorting and projection.

for the covariance, we rotate it by the 3×3 part of the view matrix W:

$$\Sigma_{view} = W \cdot \Sigma \cdot W^T$$

let w_mat = Mat3::from_mat4(view);
let w_mat_t = w_mat.transpose();

let cov3d_view = w_mat * cov3d * w_mat_t;

this is just the standard basis-change formula for a covariance matrix. the shape of the ellipsoid does not change due to this infact we are only re-expressing it in the camera's coordinate system.

1.3: projecting to 2D

now we have a 3D Gaussian in view space and we need to project it onto the 2D image plane. the projection is perspective, which means a 3D Gaussian does not project to an exact 2D Gaussian because perspective is a nonlinear transform. but we can locally linearize it using the Jacobian of the projection function, and the result is close enough.

the projection function maps a 3D point (x, y, z) in view space to pixel coordinates (u, v):

$$u = f_x \cdot \frac{x}{z_c} + c_x$$ $$v = f_y \cdot \frac{y}{z_c} + c_y$$

where $f_x, f_y$ are the focal lengths and $c_x, c_y$ are the principal point (image center). for simplicity of calcuations we can also assume, $f_x = f_y$

the Jacobian J of this projection evaluated at the splat center is:

$$J = \begin{bmatrix} \frac{f_x}{z_c} & 0 & \frac{f_x \cdot x_v}{z_c^2} \ 0 & \frac{f_y}{z_c} & \frac{f_y \cdot y_v}{z_c^2} \ 0 & 0 & 0 \end{bmatrix}$$

the structure of this matrix is very sparse, only 4 of the 9 entries are nonzero.we can just do the full $JC J^T$ with two 3×3 matrix multiplies (~54 scalar multiplies). but because we only need the top-left $2\times 2$ of $J$ Cov3D_view $J^T$, since the third row of $J$ is all zeros. and the first two rows of J each have only two nonzero entries. so instead of two full matrix multiplies, we can compute the 2D covariance with ~20 scalar multiplies by expanding the product by hand:

---

let c = &cov3d_view;
let inv_zc = 1.0 / zc;
let inv_zc2 = inv_zc * inv_zc;

let j00 = fx * inv_zc;
let j02 = fx * xv * inv_zc2;
let j11 = fy * inv_zc;
let j12 = fy * yv * inv_zc2;

// Row 0 of J * C: [j00*c00 + j02*c20, j00*c01 + j02*c21, j00*c02 + j02*c22]
let t0x = j00 * c.x_axis.x + j02 * c.z_axis.x;
let t0y = j00 * c.y_axis.x + j02 * c.z_axis.y;
let t0z = j00 * c.x_axis.z + j02 * c.z_axis.z;

// Row 1 of J * C: [j11*c10 + j12*c20, j11*c11 + j12*c21, j11*c12 + j12*c22]
let t1y = j11 * c.y_axis.y + j12 * c.z_axis.y;
let t1z = j11 * c.y_axis.z + j12 * c.z_axis.z;

// 2D cov = (J*C) * J^T, top-left 2x2:
let cov2d_00 = t0x * j00 + t0z * j02 + eps2d;
let cov2d_01 = t0y * j11 + t0z * j12;
let cov2d_11 = t1y * j11 + t1z * j12 + eps2d;

notice the eps2d on the diagonal entries. that is a small dilation (default 0.3) added for numerical stability, it just ensures the 2D covariance is strictly positive definite (not just semi-definite), which means it is always invertible.

---

Note: why the eps2d trick works

by construction, the 2D covariance $JCJ^T$ is only positive semi-definite ($A^T A$ form). but we need to invert it later (for evaluating the Gaussian at each pixel). a singular matrix is not invertible.

adding eps2d to the diagonal means adding $\lambda I$ to the matrix. for any vector x:

$$x^T \cdot (A^T A + \lambda I) \cdot x = |Ax|^2 + \lambda |x|^2 > 0$$

this is strictly positive for any nonzero x, which is the definition of positive definite, invertible, with all eigenvalues strictly positive

---

next we invert the $2\times 2$ covariance. for a $2\times 2$ matrix the inverse has a closed form:

$$\begin{bmatrix} a & b \ b & d \end{bmatrix}^{-1} = \frac{1}{ad - b^2} \begin{bmatrix} d & -b \ -b & a \end{bmatrix}$$

let det = cov2d_00 * cov2d_11 - cov2d_01 * cov2d_01;
if det <= 0.0 {
    return None;
}

let inv_det = 1.0 / det;
let cov2d_inv = Mat2::from_cols(
    Vec2::new(cov2d_11 * inv_det, -cov2d_01 * inv_det),
    Vec2::new(-cov2d_01 * inv_det, cov2d_00 * inv_det),
);

and the screen position is a standard perspective divide:

let sx = fx * xv * inv_zc + cx;
let sy = fy * yv * inv_zc + cy;

at this point, for each splat we have: screen position (sx, sy), depth zc, and the inverse 2D covariance matrix cov2d_inv. this is everything we need to evaluate the Gaussian at any pixel on screen.

step 2: computing bounding boxes

We would now have to evalute the splat at every pixel in the screen because the 2D gaussian has infinite support, it never truly reaches zero, but we can compute a bounding box that encloses the region where the gaussian has any visible effect and only evalute pixels inside this bounding box

the original 3DGS code computes the two eigenvalues $\lambda_1$, $\lambda_2$ of the 2D covariance (the variances along the two principal axes of the ellipse), takes $r = 3\sqrt{max(\lambda_1, \lambda_2)}$ (the 3-sigma rule, covering 99.7% of the Gaussian), and uses a circle of that radius as the bounding box

this is simple but wasteful. when a Gaussian is elongated (one eigenvalue much larger than the other), the bounding circle will also include a lot of empty space

We can create a tighter bounding box by observing that:

1. the extent along each axis is $k\sqrt{\Sigma_{ii}}$ where $\Sigma_{ii}$ is the diagonal entry of the 2D covariance for that axis. for an elongated ellipse, the short-axis extent is much smaller than the long-axis extent, so the bounding box is tighter.

2. the classic $3\sigma$ rule is conservative. a faint splat (low opacity) does not need $3\sigma$ because its contribution drops below the visibility threshold much sooner. the cutoff distance $k$ can be computed when our splat's opacity falls below certain threshold $\tau$

based on this each pixel gets

$$\alpha = \text{opacity} \cdot \exp(-\tfrac{1}{2} \cdot \mathbf{d}^T \Sigma^{-1} \mathbf{d})$$

we want $\alpha \geq \tau$, which rearranges to:

$$\mathbf{d}^T \Sigma^{-1} \mathbf{d} \leq 2 \ln\left(\frac{\text{opacity}}{\tau}\right) = k^2$$

for a near-opaque splat (opacity = 1), $k^2 = 2ln(255) = 11.1$, which is close to the $3\sigma$ value of 9. for a faint splat (opacity = 0.1), $k^2 = 2ln(255) = 11.1$, so the box and computation shrinks substantially.

```rust if s.opacity <= alpha_threshold { return None; } let k2 = (2.0 * (s.opacity / alpha_threshold).ln()).min(max_k2); if !(k2 > 0.0) { return None; }

let rx_f = (k2 * cov2d_00).sqrt(); let ry_f = (k2 * cov2d_11).sqrt(); if !rx_f.is_finite() || !ry_f.is_finite() || rx_f < 1.0 || ry_f < 1.0 { return None; }

let x0 = (sx - rx_f).floor() as i32