3D Gaussian Splatting

3D Gaussian splatting (3DGS) represents a scene as a large set of semi-transparent, anisotropic 3D Gaussians and renders an image by projecting them onto the image plane, sorting them by depth and alpha blending them front to back. Like a neural radiance field, it is fitted per scene to photographs with known camera poses, and it uses the same compositing as volume rendering; unlike NeRF, it has no network and no ray marching, so it renders in real time. Kerbl, Kopanas, Leimkühler and Drettakis introduced it in ACM Transactions on Graphics in 2023 (SIGGRAPH). It is an explicit, point-based neural rendering model, a point cloud whose points have extent. This article covers the model and its variants; the 2023 paper has a separate article.

Problem Addressed

The task is novel view synthesis of captured scenes. In 2023 the best quality on unbounded real scenes came from Mip-NeRF 360, which the authors report needs up to 48 hours of training; Instant-NGP and Plenoxels trained in minutes but at lower quality and, by §1’s estimate, about 10–15 frames per second (Kerbl et al., 2023; their Table 1 measures 2.8–17 fps for these methods). None rendered such scenes in real time at 1080p; the authors blame the many costly, noisy samples that continuous representations need per ray. Neural point-based renderers were fast but depended on multi-view stereo geometry and inherited its holes and spurious surfaces (§2.3). 3DGS starts from the sparse points of structure from motion and lets the optimization create, move and delete primitives.

Architecture

The model is a set of Gaussians plus a differentiable rasterizer, with no network; 3D representations counts the 59 parameters per Gaussian and compares the memory with other representations.

Anisotropic 3D Gaussians

Each primitive is an unnormalized Gaussian with mean $\boldsymbol{\mu}$ and covariance $\Sigma$ in world space,

$$G(\mathbf{x}) = \exp\!\Big(-\tfrac{1}{2} (\mathbf{x} - \boldsymbol{\mu})^\top \Sigma^{-1} (\mathbf{x} - \boldsymbol{\mu})\Big),$$

which the paper’s Equation (4) prints with $\mathbf{x}$ in place of $\mathbf{x} - \boldsymbol{\mu}$, although the text says the Gaussian is centred at $\boldsymbol{\mu}$. A covariance must be positive semi-definite, and the authors note that gradient descent on its six free entries easily leaves that set (§4). They therefore store a scale vector $\mathbf{s} \in \mathbb{R}^3$ and a quaternion $\mathbf{q}$ and build (Eq. 6)

$$\Sigma = R\,S\,S^\top R^\top, \qquad S = \operatorname{diag}(\mathbf{s}), \qquad R = R(\mathbf{q}/\lVert \mathbf{q} \rVert).$$

With $M = RS$, $\mathbf{v}^\top \Sigma \mathbf{v} = \lVert M^\top \mathbf{v} \rVert^2 \ge 0$ for every $\mathbf{v}$, so any parameter values give a valid covariance; the eigenvectors are the columns of $R$, the ellipsoid’s axes, and the eigenvalues $s_k^2$ (our derivation). The scales pass through an exponential activation and the opacity $o$ through a sigmoid (§5.1), so $\Sigma$ is in fact positive definite. The paper’s Appendix A writes the rotation in terms of a unit quaternion with the real part first (Eq. 10), the $(w, x, y, z)$ order of rotation representations, which the official code also uses.

View-dependent colour

Colour is a function of the viewing direction $\mathbf{d}$, stored per Gaussian as real spherical-harmonic coefficients $k_{\ell m}$ for each colour channel,

$$c(\mathbf{d}) = \sum_{\ell=0}^{3} \sum_{m=-\ell}^{\ell} k_{\ell m}\, Y_\ell^m(\mathbf{d}),$$

four bands (§7.1), that is, degree 3 (the official code’s default sh_degree), with 16 coefficients per channel and 48 in all. In the official code, $\mathbf{d}$ is the unit vector from the camera centre to the Gaussian’s mean, and the result is offset by 0.5 and clamped at zero.

Projection: EWA splatting

Each Gaussian is projected to a 2D Gaussian on the image plane, following Zwicker et al. (2001). With the world-to-camera pose $\mathbf{X}_c = R_c \mathbf{X}_w + \mathbf{t}$ of camera parameters, the mean goes to $\mathbf{t}_c = (t_x, t_y, t_z)^\top = R_c \boldsymbol{\mu} + \mathbf{t}$, and the pinhole map $\pi(\mathbf{t}_c) = (f_x t_x / t_z + c_x,\; f_y t_y / t_z + c_y)$ gives the 2D mean $\boldsymbol{\mu}'$. Perspective projection is not affine, so a Gaussian does not stay Gaussian; EWA splatting replaces the projection by its first-order Taylor expansion at the mean (Zwicker et al., 2001, Eqs. 17–20), under which the covariance transforms linearly:

$$\Sigma' = J\,W\,\Sigma\,W^\top J^\top, \qquad J = \frac{\partial \pi}{\partial \mathbf{t}_c} = \begin{pmatrix} f_x / t_z & 0 & -f_x t_x / t_z^2 \\ 0 & f_y / t_z & -f_y t_y / t_z^2 \end{pmatrix},$$

where $W$ is the rotation part of the viewing transformation, here $R_c$. This is the paper’s Equation (5), which forms a $3 \times 3$ matrix and drops its third row and column. The footprint is $G'(\mathbf{p}) = \exp(-\frac{1}{2}(\mathbf{p} - \boldsymbol{\mu}')^\top \Sigma'^{-1} (\mathbf{p} - \boldsymbol{\mu}'))$. Zwicker et al. note that the affine map is exact only for the ray through the centre, with errors growing with the divergence of the rays crossing the kernel; the example below measures them.

The official code adds 0.3 to the diagonal of $\Sigma'$, a screen-space dilation that keeps every splat’s standard deviation above about half a pixel; Yu et al. (2024) note that the paper does not mention it and trace aliasing artifacts to it. The code also clamps $t_x / t_z$ and $t_y / t_z$ to 1.3 times the tangent of the half field of view before computing $J$, which keeps the Jacobian bounded far outside the frustum; the paper describes only a guard band that rejects such Gaussians (§6).

Rasterization and blending

A pixel’s colour is the front-to-back blend of the splats that overlap it, sorted by depth (Eq. 3):

$$C(\mathbf{p}) = \sum_{i} \mathbf{c}_i\, \alpha_i \prod_{j<i} (1 - \alpha_j), \qquad \alpha_i = o_i\, G'_i(\mathbf{p}).$$

This is the compositing sum of volume rendering, $\hat C = \sum_i T_i \alpha_i \mathbf{c}_i$ with $T_i = \prod_{j<i}(1 - \alpha_j)$: a splat acts like a ray segment of optical depth $\sigma_i \delta_i = -\ln(1 - \alpha_i)$ (our reading). The authors make the same argument (Eqs. 1–3); their Equation (2) prints $\alpha_i$ for $\alpha_j$ in the product, in both the authors’ TOG version and arXiv v1.

The rasterizer is built for speed on the GPU (§6 and Appendix C):

1. Split the image into 16 × 16 pixel tiles; cull Gaussians whose 99% confidence ellipsoid misses the view frustum.
2. Duplicate each projected Gaussian once per tile it overlaps, with a 64-bit key whose high bits hold the tile index and low 32 bits the view-space depth.
3. Sort all keys with one GPU radix sort; each tile’s list is then a contiguous, depth-sorted range.
4. Launch one thread block per tile; each pixel walks its tile’s list front to back, accumulating colour and opacity, and stops once the pixel is saturated.

Sorting is per tile by each Gaussian’s centre depth, not per pixel, so the blend is approximate where splats overlap in depth; the authors argue this becomes negligible as splats approach pixel size. Appendix C adds numerical safeguards: an $\alpha$ below 1/255 is skipped, $\alpha$ is clamped at 0.99, and blending stops before the accumulated opacity would exceed 0.9999, that is, before $T$ drops below $10^{-4}$, the early ray termination of volume rendering.

The rasterizer is differentiable, and unlike Pulsar, the sphere rasterizer that inspired it, it lets every splat in a pixel receive gradients. The backward pass walks each tile’s list back to front and recovers intermediate transmittances from the final one by division, so it stores one number per pixel (§6). Limiting gradients to the 10 front-most splats, as in the “Limited-BW” ablation, cost the Truck scene about 11 dB of PSNR at 30,000 iterations (Table 3).

Training

The input is a set of images of a static scene with cameras calibrated by structure from motion (the paper cites Schönberger and Frahm’s COLMAP method), whose sparse point cloud initializes one Gaussian per point (§3). Each Gaussian starts isotropic, with a radius equal to the mean distance to its three nearest points (§5.1). On the synthetic Blender scenes, 100,000 random Gaussians in the scene’s bounding volume suffice (§7.2).

Each iteration renders one training view and steps Adam (Algorithm 1) on

$$\mathcal{L} = (1 - \lambda)\,\mathcal{L}_1 + \lambda\,\mathcal{L}_{\text{D-SSIM}}, \qquad \lambda = 0.2,$$

(Eq. 7), where D-SSIM is a dissimilarity derived from SSIM, the structural similarity index described in the neural rendering hub. The position learning rate decays exponentially. Optimization warms up at a quarter of the image resolution and upsamples twice, after 250 and 500 iterations. Spherical harmonics start at band 0 only and gain one band every 1,000 iterations until all four are active, because when a capture misses some viewing directions, optimizing all bands at once can produce wrong base (zero-order) colours (§7.1).

Adaptive density control

After the warm-up, every 100 iterations the model (§5.2):

• densifies each Gaussian whose average magnitude of view-space positional gradient exceeds $\tau_{\text{pos}} = 0.0002$. The authors observe that both missing geometry and Gaussians covering too much of the scene produce large gradients. A small Gaussian in an under-reconstructed region is cloned: a copy of the same size is moved along the gradient. A large one in an over-reconstructed region is split into two, with scales divided by $\phi = 1.6$ and positions sampled from the original as a probability density. Cloning adds volume; splitting conserves it.
• prunes Gaussians whose opacity has fallen below a threshold $\epsilon_\alpha$, and, periodically, those that are very large in world space or have a large footprint in the image.

Every 3,000 iterations all opacities are set close to zero. Gaussians that matter regain opacity and the rest are pruned, which moderates the growth in Gaussians that floaters near the cameras cause. The introduction reports 1–5 million Gaussians “for all scenes tested” (§1), but §7.2 gives about 200,000–500,000 for the synthetic scenes, so the range can apply only to the real scenes. The paper leaves some values to the code: in the official repository’s defaults (commit 54c035f), densification runs from iteration 500 to 15,000, $\epsilon_\alpha = 0.005$, the reset caps opacity at 0.01, and a Gaussian is split rather than cloned when its largest scale exceeds 1% of the scene extent.

Training runs for 30,000 iterations, with a 7,000-iteration configuration as a faster alternative (§7.1). In the ablation (Table 3), the full model’s average PSNR of 26.05 at 30,000 iterations fell to 23.90 without splitting, 25.23 with isotropic Gaussians, 25.35 without spherical harmonics, and 20.42 from random rather than SfM initialization.

Inference

Rendering a new view runs only the forward rasterizer. The synthetic scenes rendered at 180–300 frames per second (§7.2), measured in an interactive viewer built on SIBR (§7.1).

Performance

All figures below are the authors’ (Kerbl et al., 2023), measured on an A6000 GPU except Mip-NeRF 360 (trained on four A100s, counted as 48 GPU-hours), with every eighth image of the real scenes held out for testing. The rendering speed and memory on Mip-NeRF 360 are quoted in 3D representations; the paper article covers the rest.

• Tanks and Temples (Truck and Train, Table 1): PSNR 23.14, SSIM 0.841, LPIPS 0.183 after 30,000 iterations and 26 min 54 s of training, against 22.22, 0.759 and 0.257 for Mip-NeRF 360 trained for 48 hours, and 21.92 PSNR for Instant-NGP (Big).
• Deep Blending (Playroom and DrJohnson, Table 1): PSNR 29.41 against 29.40 for Mip-NeRF 360.
• Synthetic Blender scenes (Table 2, random initialization): mean PSNR 33.32, against 33.09 for Mip-NeRF as reported in its paper.

The text gives the training time as 35–45 minutes (§7.2), but Table 1 lists 41 min 33 s, 26 min 54 s and 36 min 2 s for the three datasets, so the Tanks and Temples time is below that range. On its own dataset Mip-NeRF 360 kept a higher PSNR (27.69 against 27.21, Table 1, a value copied from its paper), while SSIM and LPIPS favoured 3DGS.

Strengths

• Real-time rendering of unbounded real scenes, at quality comparable to the best NeRF variant of 2023.
• Training in minutes to tens of minutes on one GPU.
• Explicit and editable: Gaussians can be moved, deleted or composited, and adapt their density to the scene.
• No network at render time: a model is a list of primitives; nerfstudio’s documentation notes that splats exported as a .ply file load in a variety of web viewers.

Limitations

• Memory. Peak GPU memory when training large scenes can exceed 20 GB in the authors’ prototype, and rendering needs the full model, several hundred megabytes for large scenes, plus 30–500 MB for the rasterizer (§7.4).
• Popping. Splats are sorted by their centres, so overlapping Gaussians can swap order as the camera moves and the image changes abruptly; the guard band contributes too (§7.4). The example below shows it.
• Elongated and splotchy artifacts. Anisotropic Gaussians can produce elongated or blotchy artifacts, especially in poorly observed regions and views far from the training cameras (§7.4, Figs. 11–12).
• Aliasing. Rendering at a different sampling rate than training, by zooming or moving closer or farther, produces erosion, dilation and high-frequency artifacts (Yu et al., 2024).
• Geometry. The Gaussians are a radiance field, not a surface; the authors leave mesh reconstruction to future work (§8). Meshes from TSDF fusion of its rendered depth maps are much less accurate than those of NeuS on DTU (Huang et al., 2024; numbers under Variants).
• Inputs. Like NeRF, it assumes a static scene, consistent lighting and accurate poses.

Historical Importance

In 2023, 3DGS joined point-based graphics, with the EWA projection of elliptical Gaussians (Zwicker et al., 2001), and radiance fields fitted by differentiable compositing. Its authors argue that it shows a continuous representation is not strictly necessary for fast, high-quality radiance-field training (§8).

Variants

Each method below changes one part of the model.

Anti-aliasing: Mip-Splatting (CVPR 2024)

Yu et al. (2024) attribute the zoom artifacts to two causes: nothing limits how small a 3D Gaussian can become, and the screen-space dilation (0.3 pixel² in the code) hides degenerate Gaussians at the training resolution. A 3D smoothing filter convolves each Gaussian with an isotropic Gaussian sized by the highest sampling rate $\hat\nu_k = \max_n f_n / d_n$ among the training cameras that see it (Eqs. 6–9), so that it carries no frequencies above the Nyquist limit of its best view. A 2D Mip filter replaces the dilation with a normalized Gaussian approximating the pixel’s box filter (Eq. 10). Trained at full resolution on Blender and rendered at one-eighth resolution, 3DGS scored PSNR 17.69 and Mip-Splatting 28.67 (Table 2); trained on Mip-NeRF 360 images downsampled eightfold and rendered at 8× that resolution, 19.59 against 26.22 (Table 3). The official code now offers the 2D filter as an option (--antialiasing), off by default.

Surfaces: 2D Gaussian Splatting and SuGaR (2024)

2D Gaussian splatting (Huang et al., SIGGRAPH 2024) flattens each primitive into an oriented elliptical disc, which has a well-defined normal. Instead of the affine projection, it intersects each pixel ray with the disc’s plane, which is exact in perspective, and adds a depth-distortion loss that concentrates blending weights along the ray and a normal-consistency loss against the normals of the rendered depth. Meshes come from TSDF fusion of rendered depth. On 15 DTU scenes it reports a mean Chamfer distance of 0.80 after 30,000 iterations and 10.9 minutes, against 1.96 for TSDF-fused 3DGS depth, 1.33 for SuGaR and 0.84 for NeuS (Table 1, authors’ arXiv v3).

SuGaR (Guédon & Lepetit, CVPR 2024) keeps 3D Gaussians but adds a regularization term that pushes them to be flat, well spread and aligned with the surface: it compares a signed distance function derived from the Gaussians’ density with the one they would have if they had those properties. It then samples points on a level set of the density and meshes them by Poisson surface reconstruction rather than marching cubes. An optional stage binds Gaussians to the mesh triangles, so editing the mesh moves them.

Compactness

Lee et al. (CVPR 2024) learn a mask that removes Gaussians contributing little, replace spherical harmonics by a hash-grid colour field, and vector-quantize scale and rotation with learned codebooks. On Mip-NeRF 360 they report 48.8 MB at PSNR 27.08, and 29.1 MB at 27.03 after quantization and entropy coding, against 746 MB at 27.46 for their re-run of 3DGS (Table 1).

Dynamic scenes: 4D Gaussian splatting (CVPR 2024)

Wu et al. (2024) keep one canonical set of 3D Gaussians and learn a deformation field: features from multi-resolution HexPlane-style voxel planes over space and time feed a small MLP that predicts each Gaussian’s change of position, rotation and scale at a given time. They report 82 frames per second at 800 × 800 on an RTX 3090 for synthetic scenes.

Feed-forward Gaussians: pixelSplat (CVPR 2024)

pixelSplat (Charatan et al., 2024) predicts Gaussians from two images in one forward pass, trained across many scenes instead of optimized per scene. An epipolar transformer matches the views, and each pixel predicts a distribution over depth from which a Gaussian’s mean is sampled, with opacity set to the sampled bucket’s probability, which makes sampling differentiable and, the authors argue, avoids local minima. On RealEstate10k and ACID image pairs it outperformed light-field transformers while rendering about 2.5 orders of magnitude faster, according to the authors.

SLAM with Gaussians

Gaussians have also become a map representation for SLAM. SplaTAM (Keetha et al., CVPR 2024) maps with isotropic Gaussians from an RGB-D camera, tracking each frame by minimizing depth and colour errors over pixels that a rendered silhouette marks as already mapped. Gaussian Splatting SLAM (Matsuki et al., CVPR 2024) works from a single monocular camera, live at 3 frames per second, optimizing poses directly against the Gaussians instead of relying on offline SfM.

Implementations

The official code is the graphdeco-inria/gaussian-splatting repository (PyTorch, a CUDA rasterizer and an SIBR-based viewer); as of October 2026 the latest commit on its default branch is 54c035f (30 October 2024). Its Gaussian-Splatting License (Inria and the Max Planck Institute for Informatics) allows non-commercial research and evaluation use only; commercial use needs the licensors’ prior consent. gsplat, from the nerfstudio project, is an independent CUDA rasterization library under the Apache 2.0 licence; as of October 2026 its latest release is v1.5.3 (4 July 2025), and its README claims the official results with up to 4× less GPU memory. Nerfstudio’s Splatfacto (ns-train splatfacto) is built on gsplat; its documentation says it will drift from the original as features are added, and that it initializes from COLMAP points when available.

Projection and its linearization error

The examples use a pinhole camera with $\mathbf{X}_c = R\,\mathbf{X}_w + \mathbf{t}$ and run in sequence as one NumPy script (about 2 s). The first block checks the factorized covariance, then compares the EWA covariance with the exact perspective projection of two million samples of a Gaussian whose largest standard deviation is 0.1; the last column, the same samples through the affine map, is sampling noise alone.

import numpy as np

def quat_to_R(q):
    """Unit quaternion (w, x, y, z), scalar first, to a rotation matrix (the paper's Eq. 10)."""
    w, x, y, z = q / np.linalg.norm(q)
    return np.array([[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)]])

def covariance(log_scale, q):
    """Sigma = R S S^T R^T with S = diag(exp(log_scale)), as stored by 3DGS (Eq. 6)."""
    M = quat_to_R(q) @ np.diag(np.exp(log_scale))
    return M @ M.T

def project(mu_w, Sigma_w, R, t, fx, fy, cx, cy):
    """Mean and 2D covariance (pixels) of a 3D Gaussian under X_c = R X_w + t and a pinhole camera."""
    x, y, z = R @ mu_w + t
    J = np.array([[fx / z, 0.0, -fx * x / z**2],          # Jacobian of (fx x/z + cx, fy y/z + cy)
                  [0.0, fy / z, -fy * y / z**2]])
    return np.array([fx * x / z + cx, fy * y / z + cy]), J @ R @ Sigma_w @ R.T @ J.T, z

rng = np.random.default_rng(0)
# Any unconstrained (log_scale, q) gives a valid covariance; a raw additive update need not.
min_eig = min(np.linalg.eigvalsh(covariance(rng.normal(-2, 1, 3), rng.normal(size=4))).min() for _ in range(10000))
Sig = covariance(np.log([0.3, 0.01, 0.01]), np.array([1.0, 0, 0, 0]))
bad = 0
for _ in range(10000):
    E = rng.normal(scale=1e-3, size=(3, 3))
    bad += np.linalg.eigvalsh(Sig + (E + E.T) / 2).min() < 0
print(f"factorized: smallest eigenvalue over 10,000 random parameter sets = {min_eig:.1e}")
print(f"raw update of a thin Gaussian (step 1e-3): {bad / 100:.1f}% of steps leave the PSD cone")

# Linearization error: exact perspective projection of samples vs J W Sigma W^T J^T.
f, W_IMG = 500.0, 640                                    # 640 px wide, focal 500 px: about 65 degrees across
R_cam, t_cam = np.eye(3), np.zeros(3)
Sig = covariance(np.log([0.10, 0.06, 0.03]), np.array([0.9, 0.3, -0.2, 0.25]))
L = np.linalg.cholesky(Sig)
Z = rng.normal(size=(2_000_000, 3))
print("depth  off-axis  sigma/z   rel. error of 2D cov   mean shift (px)   sampling noise")
for depth, ang in [(0.5, 0), (1, 0), (2, 0), (4, 0), (8, 0), (2, 15), (2, 25), (2, 32)]:
    mu = np.array([depth * np.tan(np.radians(ang)), 0.0, depth])
    m2, S2, _ = project(mu, Sig, R_cam, t_cam, f, f, W_IMG / 2, W_IMG / 2)
    X = mu + Z @ L.T
    uv = np.stack([f * X[:, 0] / X[:, 2], f * X[:, 1] / X[:, 2]], 1) + W_IMG / 2
    x, y, z = mu
    J = np.array([[f / z, 0, -f * x / z**2], [0, f / z, -f * y / z**2]])
    lin = (X - mu) @ J.T + m2                            # the same samples through the affine approximation
    rel = lambda A: np.linalg.norm(A - S2) / np.linalg.norm(S2)
    print(f"{depth:5.1f}  {ang:5d} deg   {0.1 / depth:.3f}     {rel(np.cov(uv.T)):.2e}"
          f"              {np.linalg.norm(uv.mean(0) - m2):.3f}           {rel(np.cov(lin.T)):.1e}")

Output:

factorized: smallest eigenvalue over 10,000 random parameter sets = 7.5e-06
raw update of a thin Gaussian (step 1e-3): 83.1% of steps leave the PSD cone
depth  off-axis  sigma/z   rel. error of 2D cov   mean shift (px)   sampling noise
  0.5      0 deg   0.200     8.65e-02              8.035           4.6e-04
  1.0      0 deg   0.100     1.93e-02              1.944           4.6e-04
  2.0      0 deg   0.050     4.41e-03              0.485           4.6e-04
  4.0      0 deg   0.025     8.43e-04              0.123           4.6e-04
  8.0      0 deg   0.013     2.95e-04              0.032           4.6e-04
  2.0     15 deg   0.050     3.66e-03              0.395           7.3e-04
  2.0     25 deg   0.050     3.14e-03              0.344           9.6e-04
  2.0     32 deg   0.050     2.83e-03              0.317           1.1e-03

Every factorized covariance is positive definite, while 83% of small unconstrained updates push a thin Gaussian out of the positive semi-definite cone. The error of $J W \Sigma W^\top J^\top$ falls roughly with the square of $\sigma/z$: 8.7% at $\sigma/z = 0.2$, 0.44% at 0.05, and below the sampling noise at 8 m; the samples’ mean also drifts from the projected mean, by 8 pixels at 0.5 m. Moving the Gaussian off-axis within this 65° field of view changed the error little. In our reading, the linearization is accurate for splats small relative to their depth and degrades for large Gaussians close to the camera.

Tiles, sorting and the volume-rendering form

The second block renders five Gaussians into a 64 × 64 image with the official code’s dilation, clamp, cutoff and stopping rule, once over the whole image and once through 16 × 16 tiles, each listing the splats whose $3\sigma$ square overlaps it (the official rasterizer’s radius).

H = W = 64
fx = fy = 70.0; cx = cy = (W - 1) / 2
th = np.radians(10)                                     # camera at the origin, yawed 10 degrees
R_cam = np.array([[np.cos(th), 0, -np.sin(th)], [0, 1, 0], [np.sin(th), 0, np.cos(th)]])
t_cam = np.array([0.0, 0.0, 0.0])
G = [  # mean, log scales, quaternion (w, x, y, z), opacity, RGB colour
    (np.array([-0.30, 0.10, 3.0]), np.log([0.45, 0.12, 0.10]), np.array([0.92, 0.0, 0.0, 0.38]), 0.9, [0.9, 0.2, 0.1]),
    (np.array([0.25, -0.05, 3.6]), np.log([0.25, 0.40, 0.10]), np.array([0.95, 0.2, 0.0, -0.2]), 0.8, [0.1, 0.4, 0.9]),
    (np.array([0.05, 0.25, 2.5]), np.log([0.10, 0.10, 0.10]), np.array([1.0, 0.0, 0.0, 0.0]), 0.95, [0.95, 0.9, 0.2]),
    (np.array([0.70, 0.30, 4.0]), np.log([0.30, 0.06, 0.30]), np.array([0.8, 0.0, 0.6, 0.0]), 0.7, [0.2, 0.8, 0.3]),
    (np.array([-0.10, -0.40, 3.3]), np.log([0.60, 0.05, 0.05]), np.array([0.97, 0.0, 0.0, 0.24]), 0.85, [0.6, 0.3, 0.8]),
]
BG = np.array([1.0, 1.0, 1.0])
splats = []
for mu, ls, q, o, c in G:
    m2, S2, z = project(mu, covariance(ls, q), R_cam, t_cam, fx, fy, cx, cy)
    S2 = S2 + 0.3 * np.eye(2)                          # the official code's screen-space dilation (h_var = 0.3)
    splats.append((z, m2, np.linalg.inv(S2), o, np.array(c), np.linalg.eigvalsh(S2).max()))
pix = np.stack(np.meshgrid(np.arange(W), np.arange(H)), -1).reshape(-1, 2).astype(float)

def alphas(order, px=pix):
    """Per-pixel alpha of each splat, in the given order: min(0.99, o G2D), zeroed below 1/255."""
    A = []
    for k in order:
        z, m2, P, o, c, _ = splats[k]
        d = px - m2
        a = np.minimum(0.99, o * np.exp(-0.5 * np.einsum("ni,ij,nj->n", d, P, d)))
        A.append(np.where(a < 1 / 255, 0.0, a))
    return np.array(A)                                  # (K, pixels)

def blend(order, t_min=1e-4, px=pix):
    """Front-to-back alpha blending in the given order, stopping per pixel once T would drop below t_min."""
    A = alphas(order, px)
    C, T = np.zeros((len(px), 3)), np.ones(len(px))
    for a, k in zip(A, order):
        live = T * (1 - a) >= t_min                    # as in the official rasterizer: skip a splat that saturates
        C += (live * T * a)[:, None] * splats[k][4]
        T = np.where(live, T * (1 - a), T)
    return C + T[:, None] * BG, T

depth_order = list(np.argsort([s[0] for s in splats]))
img, T_final = blend(depth_order)

# Tile-based version: 16 x 16 tiles, each splat listed in the tiles its 3-sigma square overlaps.
tiles = {}
for k in depth_order:                                   # depth-sorted, so each tile list is sorted too
    z, m2, P, o, c, lam = splats[k]
    r = np.ceil(3 * np.sqrt(lam))
    for ty in range(max(0, int((m2[1] - r) // 16)), min(H // 16, int((m2[1] + r) // 16) + 1)):
        for tx in range(max(0, int((m2[0] - r) // 16)), min(W // 16, int((m2[0] + r) // 16) + 1)):
            tiles.setdefault((ty, tx), []).append(k)
img_t = np.tile(BG, (len(pix), 1))
for (ty, tx), lst in tiles.items():
    m = (pix[:, 0] // 16 == tx) & (pix[:, 1] // 16 == ty)
    img_t[m] = blend(lst, px=pix[m])[0]
print(f"{len(splats)} Gaussians, {len(tiles)} of 16 tiles touched, {sum(map(len, tiles.values()))} tile-splat pairs")
print(f"tiled vs whole-image blending: max |difference| = {np.abs(img_t - img).max():.1e}")
print(f"pixels where the splats block over 1% of the background (T < 0.99): {(T_final < 0.99).mean():.0%}")

Output:

5 Gaussians, 16 of 16 tiles touched, 53 tile-splat pairs
tiled vs whole-image blending: max |difference| = 3.4e-03
pixels where the splats block over 1% of the background (T < 0.99): 35%

Tiling evaluates 53 tile–splat pairs instead of 80 and changes the image by at most 0.0034, presumably where a splat’s tail beyond its $3\sigma$ square still has an $\alpha$ above 1/255; the official rasterizer makes the same cut.

The third block reorders the splats, recomputes the image in the volume-rendering form, tests the stopping rule, and moves a camera past two overlapping Gaussians at the same depth.

psnr = lambda a, b: 10 * np.log10(1 / np.mean((a - b) ** 2))
rev = blend(depth_order[::-1])[0]
print(f"reversed order: max |difference| {np.abs(rev - img).max():.3f}, PSNR vs sorted {psnr(rev, img):.1f} dB")
P = [psnr(blend(list(rng.permutation(len(splats))))[0], img) for _ in range(20)]
print(f"20 random orders: PSNR vs sorted from {min(P):.1f} to {max(P):.1f} dB")

# The same pixels through the volume-rendering form: sigma_i delta_i = -ln(1 - alpha_i).
A = alphas(depth_order)
cols = np.array([splats[k][4] for k in depth_order])
tau = -np.log(1 - A)                                    # optical depth of each splat at each pixel
T_vr = np.exp(-np.concatenate([np.zeros((1, len(pix))), np.cumsum(tau, 0)]))   # T_1 .. T_{K+1}
w = T_vr[:-1] * (1 - np.exp(-tau))                      # T_i alpha_i
C_vr = w.T @ cols + T_vr[-1][:, None] * BG
C_prod, _ = blend(depth_order, t_min=0.0)               # product form, no early stop
print(f"volume-rendering form vs product form: max |difference| = {np.abs(C_vr - C_prod).max():.1e}")
print(f"weights sum to 1 - T: max error {np.abs(w.sum(0) - (1 - T_vr[-1])).max():.1e}; "
      f"early stop at T < 1e-4 changes pixels by at most {np.abs(img - C_prod).max():.1e}")

# Early termination: 12 stacked splats of opacity 0.9 over one pixel.
a_stack, c_stack = np.full(12, 0.9), rng.uniform(0, 1, (12, 3))
C, T, used = np.zeros(3), 1.0, 0
for a, c in zip(a_stack, c_stack):
    if T * (1 - a) < 1e-4:                              # the official rasterizer's stopping test
        break
    C, T, used = C + T * a * c, T * (1 - a), used + 1
T_all = np.concatenate([[1.0], np.cumprod(1 - a_stack)])
C_all = (T_all[:-1] * a_stack) @ c_stack + T_all[-1] * BG
print(f"stack of 12: {used} splats blended before stopping, |error| {np.abs(C + T * BG - C_all).max():.1e} (bound T = {T:.0e})")

# Popping: two overlapping splats whose centre depths cross as the camera moves sideways.
pair = [("red", np.array([-0.05, 0.0, 3.0]), 0.9, np.array([0.9, 0.1, 0.1])),
        ("blue", np.array([0.05, 0.0, 3.0]), 0.9, np.array([0.1, 0.1, 0.9]))]
Sig_pair = covariance(np.log([0.3, 0.3, 0.3]), np.array([1.0, 0, 0, 0]))
print("camera x   in front   depth gap   centre pixel RGB")
for cam_x in [-0.4, -0.2, -0.01, 0.01, 0.2, 0.4]:
    th = np.arctan2(-cam_x, 3.0)                        # look at (0, 0, 3) from (cam_x, 0, 0)
    Rc = np.array([[np.cos(th), 0, -np.sin(th)], [0, 1, 0], [np.sin(th), 0, np.cos(th)]])
    tc = -Rc @ np.array([cam_x, 0.0, 0.0])
    out = []
    for name, mu, o, c in pair:
        m2, S2, z = project(mu, Sig_pair, Rc, tc, fx, fy, cx, cy)
        d = np.array([cx, cy]) - m2
        out.append((z, name, o * np.exp(-0.5 * d @ np.linalg.solve(S2 + 0.3 * np.eye(2), d)), c))
    out.sort(key=lambda s: s[0])                        # sort by the depth of each centre
    C, T = np.zeros(3), 1.0
    for z, name, a, c in out:
        C, T = C + T * a * c, T * (1 - a)
    C += T * BG
    print(f"{cam_x:+6.2f}     {out[0][1]:5s}     {out[1][0] - out[0][0]:.4f}      {np.round(C, 3)}")

Output:

reversed order: max |difference| 0.347, PSNR vs sorted 32.0 dB
20 random orders: PSNR vs sorted from 32.0 to 73.2 dB
volume-rendering form vs product form: max |difference| = 2.2e-16
weights sum to 1 - T: max error 2.2e-16; early stop at T < 1e-4 changes pixels by at most 0.0e+00
stack of 12: 3 splats blended before stopping, |error| 2.2e-04 (bound T = 1e-03)
camera x   in front   depth gap   centre pixel RGB
 -0.40     red       0.0132      [0.822 0.111 0.191]
 -0.20     red       0.0067      [0.822 0.111 0.191]
 -0.01     red       0.0003      [0.821 0.111 0.191]
 +0.01     blue      0.0003      [0.191 0.111 0.821]
 +0.20     blue      0.0067      [0.191 0.111 0.822]
 +0.40     blue      0.0132      [0.191 0.111 0.822]

Order matters: reversing it changes pixels by up to 0.35, and random orders score from 32 dB to 73 dB against the sorted image, depending on which splats swap. The two sums agree to $2 \times 10^{-16}$, and the weights sum to the opacity $1 - T$. No pixel of the sparse scene saturates, so stopping changes nothing there; over a stack of 0.9-opacity splats it stops after three, with an error of $2.2 \times 10^{-4}$, below the bound $T = 10^{-3}$. The last table is popping in miniature: a 2 cm camera move swaps which centre is nearer, by 0.3 mm, and the centre pixel jumps from red to blue.

Fitting Gaussians with densification

The last example fits 2D Gaussians, blended in a fixed order, to a 64 × 64 synthetic image (a disc, a rotated square and a thin bar on grey) with Adam on an L1 loss. It starts from 8 Gaussians coloured by the pixel under them, a stand-in for coloured SfM points. Every 200 steps it prunes Gaussians wider than half the image and splits the half with the largest accumulated positional gradient, sampling two children from each parent with scales divided by 1.6. Each run is compared with its final number of Gaussians trained from the start for the same 1,200 steps, over seeds 0–2 (about 30 s; two runs gave identical output).

import math, torch

torch.manual_seed(0)
torch.set_num_threads(4)
H = W = 64
yy, xx = torch.meshgrid(torch.arange(H, dtype=torch.float32), torch.arange(W, dtype=torch.float32), indexing="ij")
pix = torch.stack([xx, yy], -1).reshape(-1, 2)

# Target: a disc, a rotated square and a thin bar on the renderer's grey background.
img = torch.full((H, W, 3), 0.5)
disc = ((xx - 22) ** 2 + (yy - 24) ** 2) < 11 ** 2
u, v = (xx - 44) * math.cos(0.5) + (yy - 40) * math.sin(0.5), -(xx - 44) * math.sin(0.5) + (yy - 40) * math.cos(0.5)
square = (u.abs() < 9) & (v.abs() < 9)
bar = ((xx - yy - 8).abs() < 1.2) & (xx > 10) & (xx < 60)
img[disc], img[square], img[bar] = torch.tensor([0.95, 0.8, 0.2]), torch.tensor([0.2, 0.7, 0.3]), torch.tensor([0.9, 0.9, 0.95])
target = img.reshape(-1, 3)

def init(n):
    """n random centres, coloured by the target pixel under them (like coloured SfM points)."""
    mu = torch.rand(n, 2) * (W - 1)
    rgb = target[mu[:, 1].round().long() * W + mu[:, 0].round().long()]
    return {"mu": mu, "log_s": torch.full((n, 2), math.log(4.0)),           # centres (px), log scales
            "theta": torch.zeros(n), "logit_o": torch.zeros(n), "rgb": rgb}  # rotation, opacity logit, colour

def render(p):
    """Front-to-back blending of 2D Gaussians in index order over a grey background."""
    c, s = torch.cos(p["theta"]), torch.sin(p["theta"])
    R = torch.stack([torch.stack([c, -s], -1), torch.stack([s, c], -1)], -2)    # (n, 2, 2)
    P = (R * torch.exp(-2 * p["log_s"])[:, None, :]) @ R.transpose(1, 2)      # inverse covariance R S^-2 R^T
    d = pix[None] - p["mu"][:, None]                                           # (n, pixels, 2)
    a = torch.sigmoid(p["logit_o"])[:, None] * torch.exp(-0.5 * torch.einsum("npi,nij,npj->np", d, P, d))
    a = a.clamp(max=0.99)
    T = torch.cumprod(torch.cat([torch.ones(1, a.shape[1]), 1 - a]), 0)       # T_1 .. T_{n+1}
    return (T[:-1] * a).T @ p["rgb"] + T[-1][:, None] * 0.5

def psnr(x):
    return (-10 * torch.log10(((x - target) ** 2).mean())).item()

LR = {"mu": 0.5, "log_s": 0.05, "theta": 0.05, "logit_o": 0.05, "rgb": 0.05}   # per-parameter learning rates

def make_opt(p):
    return torch.optim.Adam([{"params": [p[k]], "lr": LR[k]} for k in p])

def fit(p, steps, every=None, rounds=0):
    """Adam on an L1 loss; every `every` steps, prune huge Gaussians and split the half with the largest gradients."""
    p = {k: v.clone().requires_grad_(True) for k, v in p.items()}
    opt, log, g = make_opt(p), [], torch.zeros(len(p["mu"]))
    for it in range(1, steps + 1):
        opt.zero_grad()
        loss = (render(p) - target).abs().mean()
        loss.backward()
        g += p["mu"].grad.norm(dim=1)                                          # accumulated screen-space gradient
        opt.step()
        if every and it % every == 0 and rounds > 0:
            log.append((len(p["mu"]), psnr(render(p))))
            with torch.no_grad():
                keep = torch.exp(p["log_s"]).max(1).values < W / 2            # prune Gaussians wider than half the image
                q, g = {k: v.detach()[keep] for k, v in p.items()}, g[keep]
                split = torch.zeros(len(g), dtype=torch.bool)
                split[g.topk(len(g) // 2).indices] = True
                reps = 1 + split.long()                                        # a split parent becomes two children
                q = {k: v.repeat_interleave(reps, 0) for k, v in q.items()}    # children keep the parent's place in the order
                child = split.repeat_interleave(reps)
                s, th = torch.exp(q["log_s"][child]), q["theta"][child]
                R = torch.stack([torch.stack([th.cos(), -th.sin()], -1), torch.stack([th.sin(), th.cos()], -1)], -2)
                q["mu"][child] += (R @ (torch.randn(len(s), 2) * s)[..., None])[..., 0]   # positions sampled from the parent
                q["log_s"][child] -= math.log(1.6)                            # scales divided by phi = 1.6
                p = {k: v.clone().requires_grad_(True) for k, v in q.items()}
            opt, g, rounds = make_opt(p), torch.zeros(len(p["mu"])), rounds - 1
    log.append((len(p["mu"]), psnr(render(p))))
    return log

print(f"PyTorch {torch.__version__}")
for seed in range(3):
    torch.manual_seed(seed)
    log = fit(init(8), steps=6 * 200, every=200, rounds=5)
    torch.manual_seed(seed)
    fixed = fit(init(log[-1][0]), steps=6 * 200)[-1][1]
    print(f"seed {seed}: " + "  ".join(f"{n}: {q:.1f}" for n, q in log) + f"   | {log[-1][0]} from the start: {fixed:.1f} dB")

Output (PyTorch 2.14.1; PSNR in dB after each stage, labelled by the number of Gaussians):

PyTorch 2.14.1
seed 0: 8: 20.4  12: 20.8  18: 25.5  25: 26.3  37: 26.8  55: 27.0   | 55 from the start: 30.6 dB
seed 1: 8: 20.2  12: 20.6  18: 25.2  27: 26.6  37: 27.5  55: 28.1   | 55 from the start: 29.8 dB
seed 2: 8: 19.0  10: 17.3  15: 17.5  22: 19.2  31: 19.4  46: 19.6   | 46 from the start: 27.6 dB

PSNR rises as Gaussians are added, from about 20 dB with 8 to 27–28 dB with 55 for seeds 0 and 1; seed 2 lost a Gaussian to pruning in the first round and stalled near 19.6 dB. In this toy, densifying never beat starting with the final number spread at random (27.6–30.6 dB). In our reading, random placement already covers a small image evenly; adaptive density control earns its keep when the initial points, like SfM points, are sparse and uneven, and when the number a scene needs is unknown.