How a computer rebuilds a real place out of soft colored dots called Gaussian Splats

A photo is a grid of numbers

To a computer a photograph is not a scene — it's a rectangular grid of tiny squares called pixels. Each pixel holds three numbers: how much red, green, and blue it shows, each from 0 (none) to 255 (full). A 1920×1080 photo is just 1920 × 1080 × 3 ≈ 6.2 million numbers. Nothing in the image is labeled; there is no "building" or "sky", only numbers.

Many algorithms don't care about color, only brightness, so they first collapse the three channels into one grayscale value per pixel — a weighted average that matches how the eye perceives lightness:

The useful information is in the changes

A patch of clear sky is thousands of nearly-identical numbers. It's useless for matching photos, because every sky pixel looks like every other — you could never tell which sky pixel in photo A is the "same" one in photo B. What is useful is where the numbers change: an edge, a corner, a speckle of texture. The change from one pixel to its neighbor is called the gradient, and almost every step in this guide is, underneath, a hunt for informative gradients.

Gradient descent: guess, measure, nudge, repeat

Both Structure from Motion (finding the cameras) and splat training (building the scene) work from the same idea: make a guess, measure how wrong it is (the "error" or "loss"), and adjust the guess so the error gets a little smaller. Repeat thousands of times and the guess gradually becomes a solution. This is gradient descent.

The picture to keep in mind: you're standing on a hillside in fog, and the height of the ground is the error. You can't see the valley, but you can feel which way the ground slopes under your feet. Step downhill a little, feel again, step again. "Rolling the error downhill" is the one engine that drives everything in this guide — bundle adjustment later in the first half, and the whole training loop in the second.

Note the word gradient is doing double duty. A moment ago it meant how much brightness changes from one pixel to the next. Here it means how much the error changes when you change one of your guessed numbers a hair. It's the same idea — "if I move this input a little, how much and which way does the output move?" — just applied to a different input. We'll also need a word for a number's share of the gradient: how much this number is responsible for the error, and which way it should move. We'll call that its blame.

Before you can rebuild a scene from photos you must know, for each photo, the exact position and orientation of the camera that took it — its pose — plus the camera's focal length (how zoomed-in it is). Get these wrong and everything downstream is built on sand. But a photo doesn't record its own pose; you only have a pile of overlapping pictures. Recovering the poses and a rough 3D point cloud purely from the pictures is Structure from Motion (SfM).

It works because the same physical thing — a corner, a sign, a rooftop — appears in many photos, and its position in each one constrains how the cameras must be arranged. The pipeline is: find distinctive spots → match them across photos → solve the two-view geometry → grow to all photos → polish everything jointly. We'll build each stage from scratch.

You need points you can recognize again in another photo taken from a different angle, distance, and rotation. The gold-standard method is SIFT (Scale-Invariant Feature Transform). It does two jobs: find good points (the detector) and describe each so it can be matched (the descriptor).

Why "scale" is the hard part, and how blur solves it

The same corner is a handful of pixels in a far photo and hundreds in a near one. If you looked for corners at a single, fixed pixel size you'd find it in one photo and miss it in the other. The trick is to look at the image at many levels of blur at once. Blurring throws away fine detail; a heavily blurred image only shows large structures. So a stack of increasingly blurred copies is really a stack of the image "seen at different sizes." SIFT searches that whole stack, so a feature is found at whatever blur level (= size) it lives.

Blurring is done with a Gaussian blur: replace each pixel with a weighted average of its neighbors, weights following the bell curve (near neighbors count more). Blur a bit, blur more, and once the image is blurry enough you can safely halve its resolution and keep going cheaply — each halving is called an octave.

// A Gaussian blur is "average each pixel with its neighbors, bell-curve weighted."
// Done in two 1-D passes (horizontal then vertical) — far cheaper than one 2-D pass.
function gaussianBlur(src, w, h, sigma) {
  const r = Math.ceil(sigma * 3);                  // how far the bell reaches
  const k = new Float32Array(2*r + 1);           // the weights
  for (let i = -r; i <= r; i++) k[i+r] = Math.exp(-(i*i) / (2*sigma*sigma));
  normalize(k);                                // weights sum to 1 (no brightening)
  // ...convolve rows with k, then columns with k...
}

Difference-of-Gaussians: turning blur into a blob detector

Now the key idea. Subtract two adjacent blur levels, pixel by pixel. In a flat region the two blurs are identical, so the difference is zero. But wherever there's a blob of roughly the size that sits between those two blur levels, the difference lights up. This subtraction, the Difference-of-Gaussians (DoG), is a cheap blob detector, and it responds most strongly at the scale matching a feature's size — exactly what we want.

A candidate feature is a pixel that is a peak in this DoG stack — brighter (or darker) than all 26 of its neighbors: the 8 around it at the same blur level, plus the 9 directly above and 9 below in the next blur levels. Being a peak across scale as well as space is what makes the detection size-independent.

// D[s] = blur[s+1] − blur[s]  (Difference-of-Gaussians at scale s)
for (let i = 0; i < ow*oh; i++) d[i] = G[s+1][i] - G[s][i];

// a candidate is a pixel bigger/smaller than all 26 neighbors in 3D (x, y, scale)
let isMax = true, isMin = true;
for (let dy = -1; dy <= 1; dy++)
  for (let dx = -1; dx <= 1; dx++)
    for (const dd of [D[s-1], D[s], D[s+1]]) {   // scale below, same, above
      const u = dd[neighbor];
      if (u >= v) isMax = false;
      if (u <= v) isMin = false;
    }
if (!isMax && !isMin) continue;                     // not a peak → discard

Refining and cleaning up the candidates

Raw peaks are only pixel-accurate and include many junk points. SIFT applies three fixes, in order:

1. Sub-pixel localization. The true peak usually lies between pixels. Fit a little 3D parabola to the DoG values around the peak and jump to its exact top (the math name for that fit is a second-order Taylor expansion). This also gives a refined contrast value at the true peak.
2. Contrast rejection. If that refined peak is too faint (below a threshold of about 0.007), it's noise, not structure — drop it.
3. Edge rejection. DoG also fires along straight edges, but an edge point slides freely along the edge, so it can't be pinned down and matched. The fix uses the local curvature. Along an edge the surface curves sharply one way and barely the other; on a corner it curves both ways. The ratio of the two curvatures is captured by the 2×2 Hessian (the matrix of second derivatives). If curvature is lopsided — the ratio exceeds ~10 — it's an edge, so discard it.

// second derivatives of the DoG at the candidate = the 2x2 Hessian H
const dxx = Dbm[ii + 1]  + Dbm[ii - 1]  - 2 * Dbm[ii];   // curvature in x
const dyy = Dbm[ii + ow] + Dbm[ii - ow] - 2 * Dbm[ii];   // curvature in y
const dxy = (Dbm[ii+ow+1] - Dbm[ii+ow-1] - Dbm[ii-ow+1] + Dbm[ii-ow-1]) / 4;
//         [ dxx  dxy ]
//   H  =  [ dxy  dyy ]     trace = dxx + dyy,   det = dxx*dyy - dxy^2
const tr = dxx + dyy;
const dt = dxx * dyy - dxy * dxy;
// an edge curves sharply one way, barely the other, so trace^2/det blows up;
// a corner curves both ways, so the ratio stays small. Keep only corners
// that are also above the contrast floor. (contrast comes from the sub-pixel fit.)
const EDGE_R = 10;
if (dt > 0 &&
    tr * tr / dt < (EDGE_R + 1) * (EDGE_R + 1) / EDGE_R &&
    Math.abs(contrast) >= PEAK_THRESH) {
  ok = true;   // survives -> becomes a keypoint
}

Giving each point an orientation (so rotation stops mattering)

The same corner photographed with a tilted camera is rotated. To match it anyway, SIFT gives every surviving point its own reference direction, then describes it relative to that direction — so the description reads the same however the photo is spun. The direction is found by making a histogram: look at all the gradients in a small circle around the point, and tally their directions into 36 bins (one per 10°), each vote weighted by the gradient's strength and its distance from the center. The tallest bin is the point's orientation.

const hist = new Float32Array(36);
// vote each nearby gradient into a direction bin, weighted by strength x nearness
for (each pixel (dx, dy) in a disc around the point) {
  const gx = (L[i + 1]  - L[i - 1])  / 2;
  const gy = (L[i + ow] - L[i - ow]) / 2;
  const mag = Math.hypot(gx, gy);
  const ang = Math.atan2(gy, gx);
  const wgt = Math.exp(-(dx*dx + dy*dy) / (2*sig*sig));   // nearer = heavier
  let bin = Math.round(ang / (2*Math.PI) * 36);
  bin = ((bin % 36) + 36) % 36;
  hist[bin] += mag * wgt;
}
// ...smooth the histogram, then take the peak(s):
const hmax = Math.max(...hist);
for (each bin b that is a local max && hist[b] >= 0.8 * hmax)
  addOrientation( parabolicPeak(b) );   // the dominant angle, plus any within 80%

The 128-number fingerprint

Finally, the descriptor. Around the point (rotated to line up with its orientation) lay a 4×4 grid of cells. In each cell, build a small 8-direction gradient histogram. That's 4 × 4 × 8 = 128 numbers — the fingerprint. Three touches make it robust:

• Trilinear spreading: each gradient contributes a little to neighboring cells and bins, so a feature landing near a cell boundary doesn't flip the fingerprint.
• Normalization: divide by the overall magnitude (so brighter/darker lighting doesn't matter), clamp any single value to 0.2 and re-normalize (so one glaring gradient can't dominate).
• RootSIFT: take the square root of each value. This quietly makes comparing two fingerprints behave better, so strong edges don't drown out subtle texture.

// 4×4 cells × 8 directions = 128. Each gradient is spread (trilinear) into
// neighboring cells/bins so the fingerprint is stable near boundaries.
acc[(vb*DW + ub)*DB + obn] += wgt * wu * wv * wo;

// normalize: L2 → clamp each to 0.2 → L2 again → RootSIFT (L1-normalize, sqrt)
for (i) acc[i] = Math.min(0.2, acc[i] * inv);      // no single spike dominates
for (i) out[i] = Math.sqrt(acc[i] / l1) * 512;    // RootSIFT, stored as bytes 0–255

Output of Stage 1: for each photo, a few thousand points, each a pixel location plus a 128-byte fingerprint that is stable across size, rotation, and lighting.

Two photos, each with a few thousand fingerprints. Which point in photo A is the same physical thing as which point in photo B? The answer for a given A-point is simply its nearest neighbor in B — the B-fingerprint with the smallest total difference. "Difference" is squared Euclidean distance: subtract the two 128-number fingerprints element by element, square, and sum.

Comparing every fingerprint in A against every one in B, for every pair of photos, is the single most expensive part of SfM — thousands of pairs, millions of 128-number comparisons each. Running it on the graphics card can make quick work of it. The GPU has thousands of tiny cores; the job is split so one core handles one A-fingerprint and streams B's fingerprints past it in blocks (a block is copied into fast on-chip "workgroup memory" that all cores in a group share, so it's read once and reused by 64 queries at a time).

// tileB is 64 B-fingerprints copied into shared on-chip memory,
// reused by all 64 query cores in this workgroup.
for (var bi = 0u; bi < tEnd; bi++) {
  var s = 0u;
  for (var k = 0u; k < 32u; k++) {          // 128 bytes packed into 32 words
    let x = da[k]; let y = tileB[bi*32u + k];
    let d0 = i32(x & 0xFFu) - i32(y & 0xFFu);  // difference of one byte
    let d1 = i32((x >> 8u) & 0xFFu) - i32((y >> 8u) & 0xFFu);
    s = s + u32(d0*d0 + d1*d1 + ...);        // squared, summed
  }
  if (s < best) { second = best; best = s; bestIdx = bBase + bi; }
  else if (s < second) { second = s; }        // track 2nd best for the ratio test
}
// best & second are returned; the CPU applies the ratio test + cross-check.

A run that would take minutes of ordinary code finishes in about a second. Output: for each pair of photos, a clean list of "this point here = that point there."

To reason about geometry we need a model of a camera. The pinhole model is enough: light from a 3D point travels in a straight line through a single point (the lens center) onto the flat sensor. A 3D point in the camera's own frame — measured from the camera's position, with z pointing straight out of the lens — p = (x, y, z), lands on the image at

Two views: the essential matrix, and why rays must meet

Take a point seen in both photos. From camera 1 it defines a ray out into space; from camera 2, another ray. Since both rays pass through the same physical point, they must intersect. Two arbitrary lines in space usually miss each other — they only meet if the two cameras are arranged just right. That "arranged just right" is captured by a single 3×3 table of numbers, the essential matrix E, which secretly encodes the rotation and the direction of travel between the two cameras.

Because each match gives one linear equation and E has 9 entries, 8 matches are enough to solve for E. (Why not 9? Multiplying all 9 entries by the same number leaves "= 0" still true, so the overall scale is free and only 8 numbers actually need pinning down.) Stack the 8+ equations into a matrix A and find the E that best satisfies "A E = 0." With noisy matches no E satisfies it exactly, so we want the one A shrinks closest to zero — found as the smallest eigenvector of A^TA. A real essential matrix must also have a specific shape (two equal non-zero "strengths" and one zero), so the raw solution is snapped to the nearest valid one.

Coping with wrong matches: RANSAC

Even after the ratio and cross-check tests, some matches are wrong (repeated textures, moving objects). A single bad match can wreck a least-squares solve — the standard "best fit" that minimizes the sum of squared misses, and which treats every input as trustworthy, so one liar drags the whole answer toward itself. The defense is RANSAC: instead of trusting all matches, repeatedly pick a random minimal set (8), solve E from just those, and count how many of all matches agree with it. Keep the E with the biggest agreement, and treat the disagreers as outliers to discard. "Agreement" is measured by the Sampson distance — essentially how far a match is from satisfying x₂^TE x₁ = 0. RANSAC even adjusts how many random tries it needs based on how clean the data looks.

// Each match adds one row (its epipolar equation) to the 9×9 system AᵀA.
for (const idx of indices) {
  a[0]=x2*x1; a[1]=x2*y1; a[2]=x2; a[3]=y2*x1; ...          // the 9 coefficients
  for (i) for (j) AtA[i*9+j] += a[i]*a[j];
}
const { vecs } = jacobiEigen(AtA, 9);
return enforceEssential(matCol(vecs, 9, 8));    // smallest eigenvector → E → snap to valid

// RANSAC: many random 8-match samples; keep the E with the most inliers,
// where an inlier's Sampson distance is below a pixel-scaled threshold.
for (let it = 0; it < iters; it++) {
  const est = estimateE(x1s, x2s, random8());
  let count = 0;
  for (i) if (sampsonError2(est.E, x1s[i], x2s[i]) < th2) count++;
  if (count > bestCount) { best = est; bestCount = count;
    iters = adaptiveIterationCount(count / n); }        // stop early if data is clean
}

From E the relative pose is recovered by decomposition: the math yields four possible (R, t) combinations, of which exactly one puts the triangulated points in front of both cameras (you can't photograph what's behind you). That "in front of both" check — cheirality — picks the true pose.

With the two cameras arranged, each matched point's two rays cross at one place — its 3D location. This is triangulation, the same depth trick two eyes use. Each camera contributes the statement "the 3D point X, projected through me, lands where I see it" — which is two straight-line (linear) equations in X's three unknowns. Two cameras give four equations for three unknowns; the best-fit X is, again, the smallest-eigenvector solution of the stacked system (the Direct Linear Transform).

export function triangulateN(cams, obs) {
  for (let k = 0; k < cams.length; k++) {
    const { R, t } = cams[k]; const [u, v] = obs[k];
    // two rows per view: (u·row3 − row1)·X = 0 ,  (v·row3 − row2)·X = 0
    row[0] = u*R[6]-R[0]; row[1] = u*R[7]-R[1]; row[2] = u*R[8]-R[2]; row[3] = u*t[2]-t[0]; addRow();
    row[0] = v*R[6]-R[3]; row[1] = v*R[7]-R[4]; row[2] = v*R[8]-R[5]; row[3] = v*t[2]-t[1]; addRow();
  }
  return smallestEigenvector(AtA);   // X as 4 numbers (homogeneous): divide the first 3 by the 4th for x,y,z
}

A guard rejects points with too little parallax — if the two rays are nearly parallel (the cameras barely moved relative to the point's distance), the crossing is wildly uncertain, so that point is dropped.

Two photos give a seed reconstruction. To fold in the rest, the pieces are stitched together with tracks and grown one camera at a time.

1. Build tracks. The same point may be matched across many photo pairs: A↔B, B↔C, C↔D. A "union-find" data structure chains these into a single track — "these observations, in these five photos, are all the same physical point." Tracks are what tie the whole reconstruction together.
2. Pick the starting pair. Choose two photos that share many tracks and have healthy parallax (a wide enough baseline for reliable triangulation), and seed the reconstruction from them.
3. Register the next camera (PnP). A new photo already sees tracks whose 3D positions are known. Given several "this 3D point appears at this pixel" pairs, solving for that one camera's pose is the Perspective-n-Point problem — a direct linear solve for a rough pose, wrapped in RANSAC (against wrong matches) and then refined by gradient descent to minimize its reprojection error.
4. Triangulate more, repeat. The newly-posed camera lets fresh tracks be triangulated, which lets the next camera register, and so on until every photo is placed.

// gather this image's features that hit already-triangulated tracks:
//   obj = known 3D points,  im = where they appear in the new photo
const pnp = pnpRansac(obj, im, threshold, rng);   // robust rough pose
if (pnp) poses[img] = refinePose(pnp.R, pnp.t, obj, im, 10); // polish by gradient descent (Levenberg–Marquardt, explained in Stage 5)
// then triangulate every fresh track this camera newly sees.

This is the first place the "roll a number downhill" engine appears in full. Define one error: for every sighting (a 3D point seen in a photo), how many pixels away does the point project from where it's actually observed? Sum that over all sightings. Then adjust all camera poses, all 3D points, and the shared focal length and lens distortion together, a little at a time, to shrink that sum. This joint polish is bundle adjustment.

With thousands of cameras and hundreds of thousands of points, this is solved not by brute force but by two standard tools:

• Levenberg–Marquardt — a careful gradient-descent step: near the answer it takes big confident steps, far away it takes small cautious ones, adjusting its "damping" each round based on whether the last step helped.
• The Schur complement — a bookkeeping trick. Because each 3D point only touches the cameras that see it, the equations have a special block structure. The huge number of point-unknowns can be eliminated in bulk, leaving a much smaller system in just the camera-unknowns to solve directly, then the points are recovered by back-substitution.

// project a point through a camera, with radial lens distortion
const xp = x/z, yp = y/z;                       // perspective divide → normalized
const r2 = xp*xp + yp*yp;
const D  = 1 + k1*r2 + k2*r2*r2;                // lens distortion
predicted = { u: f*xp*D + cx, v: f*yp*D + cy };

// each round: measure how every pixel miss responds to every unknown (the
// "normal equations", JᵀJ Δ = −Jᵀr — the standard least-squares step), Huber-weight
// big residuals, damp (Levenberg–Marquardt), Schur-eliminate the point block,
// solve for the camera step, apply it, and keep it only if the total error fell.

SfM handed us posed cameras and a sparse point cloud. Now we reconstruct the look of the scene — not as surfaces, but as a cloud of soft, semi-transparent, colored blobs called Gaussian splats. Enough of them, overlapping, melt into what reads as solid geometry. The reason blobs are used instead of triangles is decisive: every property of a blob is a plain number you can nudge, and nudging it changes the rendered image smoothly. Smoothness is exactly what gradient descent needs.

To score a guess we must draw the current blob cloud through one of the known cameras. This is the forward render: one projection step that flattens every blob onto the screen, followed by five tiling stages that make painting them fast. First the per-blob geometry, then the tiling.

Flattening a 3D blob into a 2D oval (EWA projection)

EWA stands for Elliptical Weighted Average — the name of the 2001 splatting technique whose projection step this is. The name matters less than the mechanic below.

Move the blob into the camera's frame (rotate and translate by the camera pose), then flatten to the image. A blob doesn't flatten to a circle — a tilted 3D ellipsoid becomes a tilted 2D ellipse, and it's stretched by perspective (nearer parts loom larger). Capturing that stretch exactly is impossible with a straight projection, so we use the local linear approximation of the projection at the blob's center — its Jacobian J, the 2×3 table of "how a nudge in each 3D direction moves the 2D image point." Push the 3D covariance through the camera rotation W and this Jacobian and you get the 2D covariance of the flattened oval.

Inverting that 2×2 gives the conic (A, B, C) — three numbers that say how fast the oval fades as you step d pixels from its center. (A matrix inverse is the "undo" of a matrix, like 1/x for a number; here it flips "how spread out the oval is" into "how quickly it fades.") A blob's contribution at a pixel is its opacity times that falloff. (A small "compensation" factor, from a refinement called Mip-Splatting, keeps very tiny blobs from adding spurious energy.)

let f = cam.proj.x; let iz = 1.0 / g.pc.z;
// rows of T = J·W, folded together
g.t0 = (f*iz)*cam.R0.xyz + (-f*g.pc.x*iz*iz)*cam.R2.xyz;
g.t1 = (f*iz)*cam.R1.xyz + (-f*g.pc.y*iz*iz)*cam.R2.xyz;
// V = T·Σ·Tᵀ  → the 2×2 screen covariance (va, vb, vc)
g.va = dot(g.t0, Sigma·t0); g.vb = dot(g.t1, Sigma·t0); g.vc = dot(g.t1, Sigma·t1);
// invert V → conic (A,B,C); opacity via sigmoid; cull if it can't reach a pixel
proj[b+3] = vc*inv; proj[b+4] = -vb*inv; proj[b+5] = va*inv;   // A, B, C
let opa = 1.0 / (1.0 + exp(-logitOpacity));

Why tiling, and the five stages

A photo has millions of pixels; the cloud has millions of blobs. Testing every pixel against every blob is hopeless. Instead the screen is cut into 16×16-pixel tiles, and each pixel only ever considers the blobs that actually land on its tile — already sorted near-to-far. Building those per-tile lists is a classic parallel pattern:

1. Count. Each blob figures out which tiles its oval covers and adds 1 to each of those tiles' counters (an atomic add, since many blobs hit the same tile at once).
2. Scan. A running total (prefix sum) of the tile counts reserves an exact, tightly-packed slot range for every tile in one big list — no wasted space, no fixed per-tile cap.
3. Scatter. Each blob writes its ID (paired with its depth) into every tile it covers, using each tile's reserved range.
4. Sort. Within each tile, sort the blobs by depth, near to far. Small tiles are sorted entirely in fast on-chip memory with a bitonic sort (a sorting network that's efficient for parallel hardware); big ones spill to a global sort.
5. Composite. Each pixel walks its tile's sorted blobs front-to-back, blending — the actual painting.

// how far from its center can this blob still tint a pixel? (in pixels)
let ad = va + 0.3;  let cd = vc + 0.3;              // Mip-dilated 2D covariance
let detVd = ad*cd - vb*vb;
let mid = 0.5*(ad + cd);
let disc = sqrt(max(mid*mid - detVd, 0.0));      // eigenvalues = mid ± disc  (the oval's long and short radii, squared)
let eMax = min(E_CUT, log(max(opacity*comp / A_MIN, 1.0001)));
let rad  = sqrt(2.0*eMax) * sqrt(mid + disc);     // σ-extent along the long axis
// then file this blob into every tile whose square meets [center ± rad].

const GS_TILE     = 16;        // tile is 16×16 px
const E_CUT       = 4.5;       // drop a pixel past this falloff (~3σ)
const A_MIN       = 0.0039;    // alpha floor; below this a blob is skipped
const SHARED_SORT = 2048;      // segments ≤ this sort in on-chip memory
const ENTRIES_CAP = 12000000;  // global (blob,tile) budget

// SCAN — exclusive prefix sum of the per-tile counts → each tile's start offset.
// (Real code does this as a parallel two-level scan in one workgroup; the result is:)
var acc = 0;
for (t in 0 .. numTiles) {
  tileStart[t] = acc;                    // tile t's segment begins here
  acc += paddedSize(tileCount[t]);       // pad big segments up to a power of two
}
tileStart[numTiles] = acc;               // sentinel closes the last segment

// BITONIC SORT — order one tile's segment [s, e) by key ascending (one workgroup/tile).
// A branch-free sorting network: build bitonic runs, then repeatedly merge them.
for (k = 2; k <= segCap; k <<= 1) {        // run length being merged
  for (j = k >> 1; j > 0; j >>= 1) {       // compare distance
    for each slot i in the segment (parallel) {
      let l = i ^ j;                       // its partner
      if (l > i) {
        let asc = (i & k) == 0;            // this block sorts ascending?
        let hi  = key[i] > key[l] || (key[i]==key[l] && id[i] > id[l]);
        if (hi == asc) { swap(i, l); }     // compare-exchange
      }
    }
  }
}
// padding slots hold key = 0xFFFFFFFF so they sink to the back and the composite stops at them.

The painting itself: front-to-back alpha compositing

Walk the sorted blobs from nearest to farthest, tracking one number, the transmittance T — "how much light is still getting through." It starts at 1 (nothing blocking) and drops as each layer covers more. Each blob adds its color scaled by how solid it is here and by how much light is still available. Once T is nearly zero the rest is hidden and the loop stops early.

var T = 1.0; var Crgb = vec3f(0.0);
for (var k = segStart; k < segEnd; k++) {          // this tile's blobs, near → far
  let d = px - vec2f(proj[b], proj[b+1]);           // pixel minus blob center
  let e = 0.5*(A*d.x*d.x + C*d.y*d.y) + B*d.x*d.y;   // oval falloff exponent
  let alpha = min(0.99, opacity * comp * exp(-e));
  if (alpha < A_MIN) { continue; }
  Crgb += T * alpha * color;                        // blend this layer in
  T *= 1.0 - alpha;                                 // less light gets through now
  if (T < 1e-4) { break; }                         // fully hidden → stop
}

Compare the rendered image to the real photo, pixel by pixel, and reduce the whole thing to one number. Per pixel, the "difference" r is the rendered color minus the photo's color. You could square it (the classic mean-squared error), but squaring punishes a few freak pixels — a moving car, glare, a blur — so hard that the cloud would contort to chase them. The trainer uses the Charbonnier loss instead: it acts like the plain absolute difference for ordinary pixels but grows only gently for huge outliers.

let err  = (gain * C + bias) - target;          // with per-photo exposure correction
const DELTA = 0.03;
let root = sqrt(err*err + DELTA*DELTA);
let loss = (root.x + root.y + root.z) - 3.0*DELTA;
let gC   = gain * (err / root);                  // slope handed to backprop

We have one loss number and its slope at each pixel. But you can't nudge a blob with a pixel's slope — you need, for every one of a blob's 14 stored numbers, "if I changed you a hair, how would the loss move?" Getting that for millions of blobs sounds impossible, but there's a mechanical answer: run the render backward. This is reverse-mode automatic differentiation, the same idea that trains neural networks, and it's just the chain rule applied with care.

The chain rule, in one breath

If A feeds into B feeds into the loss, then "how the loss depends on A" = "how the loss depends on B" × "how B depends on A." So if you know the sensitivity at the output and the small local relationship at each step, you can multiply your way backward from the loss to every input. The forward render is a chain: blob numbers → 3D covariance → 2D oval → per-pixel contribution → blended color → loss. Backprop walks that chain in reverse, carrying a "sensitivity" and multiplying by each step's local relationship.

The clever bit: a single backward sweep through the blend

The compositing has an awkward property: a blob affects a pixel two ways — it adds its own color, and it hides some of everything behind it. To get a blob's blame you'd naively need to know "everything behind it," which changes at every blob. The trick is to sweep the sorted list back-to-front, accumulating a running sum S of "everything already behind this point." Then each blob's blame comes out from a short formula — no search, no second pass.

From that one number the chain rule cascades further: onto the blob's screen position and oval shape (the conic), then — because the conic came from the 3D covariance which came from scale and quaternion — back onto its raw stored numbers. That last hop, conic → 2D covariance → 3D covariance → scale & quaternion, is a stack of matrix-calculus derivatives; the code calls it the "chain" pass.

// back-to-front through the tile; S = "everything behind this blob so far"
var S = bg * T;
for (var kk = last; kk > segStart; kk--) {
  let Tb = Ta / (1.0 - alpha);                    // light in front of this blob
  let gcolor = gC * (alpha * Tb);                 // blame for its color
  var galpha = dot(gC, c*Tb - S/(1.0 - alpha));   // blame for its solidity
  // → blame for screen position & conic (A,B,C), accumulated with atomic adds
  S += c * alpha * Tb; Ta = Tb;
}
// the "chain" pass then flows conic → 2D cov → 3D cov → scale & quaternion,
// producing the gradient on each blob's actual stored 16 numbers.

One WebGPU reality. The trainer fuses the forward paint and the backward sweep into one shader pass, so each tile's blob list is touched once. (Thousands of pixels blame the same blob at the same instant; how those blames are summed safely — the atomAdd calls and the constant 16384 you'll see below — is covered at the end of this section.)

From a pixel's blame to the oval's shape and place

The sweep above hands us, for each blob at each pixel, two blames: how its color should change and how its solidity (its alpha) should change. Color is easy — it flows straight to the blob's RGB. Alpha is the interesting one, because alpha was built from several things and the blame has to be split among them. Recall how a pixel's alpha came to be: the blob has a solid opacity, and that opacity is scaled down by how far the pixel sits from the oval's center — far pixels get almost none. "How far, oval-shaped" is one number e, computed from the three conic numbers (A, B, C) that define the oval and the pixel's offset (dx, dy) from the center:

let ga = galpha * araw;                          // through exp: ∂α/∂e = −α
let gmean = ga * vec2f(cA*d.x + cB*d.y, cB*d.x + cC*d.y);
atomAdd (0u, gmean.x); atomAdd(1u, gmean.y);      // → oval center (screen)
atomAddC(2u, -ga * 0.5 * d.x*d.x);                // → conic A  (px²-scaled, see below)
atomAddC(3u, -ga * d.x*d.y);                       // → conic B
atomAddC(4u, -ga * 0.5 * d.y*d.y);                // → conic C
atomAdd (5u, galpha * opa * G);                    // → comp
atomAdd (6u, ga * (1.0 - opa));                    // → opacity logit
atomAdd (7u, gcolor.r); atomAdd(8u, gcolor.g); atomAdd(9u, gcolor.b);

The covariance chain: oval → 3D ellipsoid → size & tilt

Now the hard hop. A, B, C describe the flat 2D oval, but a blob doesn't store an oval — it stores a 3D size and tilt. In the forward render those became a 3D covariance Σ, which the camera flattened into a 2D covariance V, which was inverted to get the conic. To assign blame we walk that same chain in reverse, one link at a time. Each link is a small matrix-calculus rule; you don't need to derive them, only to see that each is mechanical.

// 1) conic → V:  dV = −Vinv · G · Vinv   (kA,kB,kC = the conic = Vinv)
let g11 = gp[2]; let g12 = 0.5*gp[3]; let g22 = gp[4];   // symmetric conic-grad
// P = Vinv·G·Vinv, then dV = −P (off-diagonals summed)
var gva = -p00; var gvb = -(p01 + p10); var gvc = -p11;

// 2) V → Σ:  dΣ = Tᵀ·H·T   (rows of T are g.t0, g.t1 = J·W)
let h11 = gva; let h12 = 0.5*gvb; let h22 = gvc;
let u0 = h11*g.t0 + h12*g.t1;  let u1 = h12*g.t0 + h22*g.t1;
let dS00 = g.t0.x*u0.x + g.t1.x*u1.x;  // …the six unique entries of dΣ…
//    dT = 2·H·T·Σ  → feeds the position/depth gradient (next section)

// 3) Σ → M → size & tilt:  dM = 2·dΣ·M   (m0,m1,m2 = rows of M = R·diag(s))
let dM0 = 2.0*(dS00*m0 + dS01*m1 + dS02*m2);   // (dM1, dM2 likewise)
let dsv = vec3f(dM0.x*g.r0.x + dM1.x*g.r1.x + dM2.x*g.r2.x, /* y */, /* z */);
gradF[b+3u] = dsv.x*g.s.x + AREG*(ls.x - meanLogScale);   // log-scale grad (+aniso reg)
// tilt: dR = dM·diag(s) → gq = (∂R/∂q · dR) → project off the norm:
let gqr = (gq - g.q*dot(g.q, gq)) / ql;      // grad of the normalized quaternion
gradF[b+6u]=gqr.x; gradF[b+7u]=gqr.y; gradF[b+8u]=gqr.z; gradF[b+9u]=gqr.w;

let dR0 = dM0 * g.s;  let dR1 = dM1 * g.s;  let dR2 = dM2 * g.s;   // dR = dM·diag(s)
let gw = 2.0*(-qz*dR0.y + qy*dR0.z + qz*dR1.x - qx*dR1.z - qy*dR2.x + qx*dR2.y);
let gx = 2.0*( qy*dR0.y + qz*dR0.z + qy*dR1.x - 2*qx*dR1.y - qw*dR1.z + qz*dR2.x + qw*dR2.y - 2*qx*dR2.z);
let gy = 2.0*(-2*qy*dR0.x + qx*dR0.y + qw*dR0.z + qx*dR1.x + qz*dR1.z - qw*dR2.x + qz*dR2.y - 2*qy*dR2.z);
let gz = 2.0*(-2*qz*dR0.x - qw*dR0.y + qx*dR0.z + qw*dR1.x - 2*qz*dR1.y + qy*dR1.z + qx*dR2.x + qy*dR2.y);
let gq  = vec4f(gw, gx, gy, gz);
let gqr = (gq - g.q*dot(g.q, gq)) / ql;         // project off the length change; ql = |q_raw|
gradF[b+6u]=gqr.x; gradF[b+7u]=gqr.y; gradF[b+8u]=gqr.z; gradF[b+9u]=gqr.w;

Where the position gradient comes from

A blob's position collects blame from two directions. The obvious one: moving the blob slides its oval's center across the screen. The subtle one: moving the blob — especially nearer or farther — changes how the projection stretches its oval, because the flatten-map T depends on depth. That second path is why depth has such a strong pull: the projection carries 1/z and 1/z² terms, so z shows up with a steep slope.

Making the parallel sums exact: fixed-point, and the conic trick

Thousands of pixels blame the same blob at the same instant. Adding all those blames safely means atomic adds — and there's a bonus: integers add in any order with no rounding drift, so the total is bit-for-bit identical however the GPU threads happen to interleave. So each floating-point blame is multiplied by 16384, rounded to an integer, added atomically, and divided back at the end. One wrinkle: the conic blames grow with distance² and could overflow the integer for a big oval, so each conic blame is pre-divided by (1 + the oval's largest radius) in the paint pass, and multiplied back by the identical factor in the chain pass — a scaling that cancels exactly and keeps every added value near 1.

const FIXED = 16384.0;
fn atomAdd(idx: u32, v: f32) {
  atomicAdd(&gradP[idx], i32(clamp(v*FIXED, -1.0e9, 1.0e9)));   // float → int, exact sum
}
// paint pass scales conic blames by cnorm = 1/(1+λmax) of the dilated 2D cov;
// chain pass recomputes the SAME λmax and multiplies it back — cancels exactly:
gp[2] *= (1.0 + lmax); gp[3] *= (1.0 + lmax); gp[4] *= (1.0 + lmax);
// then decode every slot: value = f32(atomicLoad(gradP[k])) / FIXED, and zero it for next step.

Proving the blame is right — the gradient check

Every formula above was derived by hand, and hand-derived gradients are famously easy to get subtly wrong — a missing ½, a dropped term, a sign. There's a brute-force way to check any one of them without trusting the algebra at all: nudge the number and watch the loss. Bump one parameter up by a hair h and re-measure the loss, bump it down by h and re-measure; the loss's true slope is just the difference over 2h. If that measured slope matches the analytic blame, the derivation is right.

Now every one of a blob's numbers has a gradient — the direction that lowers the loss. The dumbest update is "subtract a little of the gradient" (plain gradient descent). It works, but badly: some numbers (position) need tiny careful steps, others (opacity) can move fast; and a number whose gradient keeps flip-flopping should slow down while one heading steadily should speed up. The optimizer that handles all this is Adam, and it keeps two running averages per number.

var m = b1*mBuf[j] + (1.0-b1)*g;               // momentum
var v = b2*vBuf[j] + (1.0-b2)*g*g;             // gradient-size memory
mBuf[j] = m; vBuf[j] = v;
let mh = m / (1.0 - pow(b1, t));               // startup bias correction
let vh = v / (1.0 - pow(b2, t));
params[j] = params[j] - lr * mh / (sqrt(vh) + eps);
// opacity also gets a gentle pull toward semi-transparent (a "regularizer":
// an extra small term in the loss that discourages a bad habit, here hoarding opacity),
// and scales/logits are clamped so they can't run off to absurd values.

Then the loop returns to step 1 with a different photo. Render → loss → backprop → Adam, about a thousand times a second, for tens of thousands of steps.

Densification: cloning the busy blobs, recycling the dead

Nudging a fixed set of blobs can't add detail the set can't represent — a fence rail, gravel, foliage — while elsewhere blobs fade to invisible and waste compute. So every couple of thousand steps the trainer reshapes the population — a densification pass. (Papers call this flavor "MCMC-style" after a random-sampling technique it borrows from; the mechanic is what matters.)

1. Find the dead and the donors. Blobs whose opacity has fallen below ~0.02 are "dead" (contributing nothing); blobs above ~0.4 are healthy "donors."
2. Relocate the dead. Rather than delete them, teleport dead blobs next to a donor (biased toward the largest donors, which usually carry the remaining blur), so the count stays roughly fixed but effort moves where it's needed.
3. Split the big ones. A large donor that keeps getting strong "not enough detail" gradients is split: both halves shrink (÷1.6) so the footprint is preserved, and they drift apart to specialize. A freshly reshaped blob has its Adam momentum reset so it can move cleanly.

Growth is front-loaded (generous early, then it stops partway through and the rest is pure polish) and hard-capped so the count can't explode.

for (let i = 0; i < n; i++) {
  const o = sigmoid(params[i*16 + 13]);          // this blob's opacity
  if (o < 0.02) dead.push(i);                    // invisible → recycle it
  else if (o > 0.4) donors.push(i);             // healthy → can spawn from it
}
// prefer splitting the biggest donors; both halves shrink so the area is kept
if (doSplit) { const shrink = Math.log(1/1.6);
  for (let k = 3; k <= 5; k++) { params[bi+k] += shrink; params[bd+k] += shrink; }
  resetAdamMoments(bd);                          // donor changed shape → forget its momentum
}

Spherical harmonics: color that changes with viewing angle

A flat-colored blob suits a matte wall, but a metal roof or wet street looks different depending on where you stand. So each blob gets a few extra color numbers that let its color shift smoothly with the viewing direction. The smooth "patterns over directions" they multiply are spherical harmonics — a fixed library of ways a value can vary over a sphere (a constant, then "brighter from above," "brighter from the left," and progressively wigglier patterns). Degree 0 is flat color; each higher degree adds finer directional detail.

Crucially this shine is unlocked gradually — one degree switched on roughly every 1000 steps — so the cloud is forced to learn what's there before how it glints, and can't cheat a hard photo by mysteriously recoloring a blob.

let v = normalize(worldPos - cameraPos);        // viewing direction to this blob
var Y = shBasis(v);                              // the fixed directional patterns
let aK = shActiveK();                            // how many are unlocked yet (ramps up)
for (var k = 0u; k < aK; k++) col += sh[..k..] * Y[k];
col = max(col, vec3f(0.0));                       // color can't go negative

So what are those "fixed patterns," concretely? Each Y_k is a simple polynomial of the viewing direction v = (x, y, z) (a unit arrow from the blob to the camera). Degree 0 is a constant — the same from every angle. Degree 1 is literally x, y, z — "more toward the right / up / front." Degree 2 and 3 are products like xy, xz, and 2z²−x²−y² — progressively wigglier ways to vary over a sphere. Each comes with a fixed numeric multiplier.

// the fixed patterns, evaluated for this view direction v=(x,y,z)
Y[0] = -0.4886025 * y;  Y[1] = 0.4886025 * z;  Y[2] = -0.4886025 * x;   // degree 1
Y[3] = 1.0925484 * x*y; Y[4] = -1.0925484 * y*z;                        // degree 2 …
Y[5] = 0.3153916 * (2.0*z*z - x*x - y*y);  Y[6] = -1.0925484 * x*z;
Y[7] = 0.5462742 * (x*x - y*y);                                         // … (degree 3 = Y[8..14])
// backward: each weight's blame is just dL/dcolor · the pattern's value
for (var k = 0u; k < activeK; k++) {
  shGrad[sb + k]          = dRGB.x * Y[k];    // red channel
  shGrad[sb + SHK + k]    = dRGB.y * Y[k];    // green
  shGrad[sb + 2u*SHK + k] = dRGB.z * Y[k];    // blue
}

The backward pass is refreshingly simple: a weight that multiplies a strongly-lit pattern in this view collects proportionally more blame — ∂L/∂sh_k = dRGB · Y_k. (This particular trainer also lets a sliver of blame flow through the view direction back onto the blob's position, which most implementations skip.) And "unlocked gradually" has a precise meaning: only the first activeK = (activeDeg+1)²−1 weights are summed, with activeDeg stepping up by one about every 1000 iterations — inactive patterns contribute nothing and receive no gradient until their band switches on.

Forgiving the camera and the light

Two last touches for real footage. SfM's poses are excellent but not perfect, so the trainer may nudge each camera's pose by a whisker as it goes (backprop already produces a gradient for the camera, since the render depends on it). And a camera's auto-exposure drifts frame to frame, so the trainer learns a small per-photo brightness gain and bias — the ones you saw in the loss — rather than twisting the blobs to satisfy inconsistent exposures.

// exposure: err = gain*C + bias - photo;  eg = err/root (the Charbonnier slope)
camAdd(camBase + 6u, dot(eg, C) * gain);      // ∂L/∂(log gain)
camAdd(camBase + 7u, eg.x + eg.y + eg.z);     // ∂L/∂bias
// pose: dpc is the cam-space position blame already computed for the blob
let torque = cross(g.pc - cam.t, dpc);
camAdd(camBase + 0u, torque.x); camAdd(camBase + 1u, torque.y); camAdd(camBase + 2u, torque.z);
camAdd(camBase + 3u, dpc.x); camAdd(camBase + 4u, dpc.y); camAdd(camBase + 5u, dpc.z);   // ∂L/∂t

Gains taper off, so the trainer keeps one photo secret — a holdout it never learns from, only scores against. While the holdout score keeps rising the cloud is genuinely learning the scene rather than memorizing the shown photos; when it flattens, training stops on its own. Because that photo was never used, its PSNR is an honest measure of quality on views the model has never seen.

What you keep is the blob cloud — every blob's final position, scale, quaternion, color (with its spherical-harmonic extras), and opacity. Fly a new camera through it in real time, from angles no photo ever captured, and it holds together; or write it to a standard .ply file that any Gaussian-splat viewer can open. Footage in, a place out.

Two solvers, one engine. Bundle adjustment nudges cameras and points to match the photos; training nudges blobs to match the photos. Both just keep rolling "how wrong" downhill.

Putting every stage together, a from-scratch implementation is:

1. Front end (SfM). SIFT features → GPU match → essential-matrix + RANSAC per pair → tracks → incremental PnP → bundle adjustment. Out: poses, focal, distortion, sparse cloud.
2. Initialize. Straighten photos with the distortion; seed 16-float blobs on the sparse points, colored from the photos; allocate gradient and Adam-moment buffers.
3. Per step. Pick a photo; project blobs (EWA); tile via count → scan → scatter → sort; composite front-to-back; Charbonnier loss; fused backward sweep + chain pass with fixed-point atomic gradients; Adam update.
4. Periodically. Every ~2500 steps, relocate dead and split big donors; every ~1000, unlock one more spherical-harmonic degree.
5. Stop & export. Watch the holdout PSNR; when it plateaus, stop and write the cloud to .ply.

Every arrow above is one of the sections you just read. There is no hidden step and no learned black box — just projection, sorting, blending, a robust error, the chain rule, and Adam, repeated until the picture stops being wrong.

And it isn't hypothetical: a working single-file trainer was built from nothing but this page — the second half, running live in your browser. It skips the SfM front end (a synthetic scene rendered from known cameras stands in for the posed photos) and trains real 16-float blobs with the exact loop above: EWA projection, tile sort, front-to-back compositing, Charbonnier loss, backprop, Adam.

This sample can generate gaussian splats from a video using the described algorithms: https://taosengine.com/samples/splat_video.html