Task 3 of 5

Lucas–Kanade: Buy a Second Equation

You cannot get more information out of one pixel, so buy it from the neighbours. Lucas–Kanade assumes that everything inside a small window moves together. A 5×5 window gives 25 copies of Ix·u + Iy·v + It = 0 sharing one unknown (u, v) — 25 equations, 2 unknowns, over-determined instead of under-determined. Least squares turns that into a 2×2 system:

| Sxx  Sxy |   | u |        | Sxt |
|          | · |   |  =  −  |     |
| Sxy  Syy |   | v |        | Syt |

where Sxx = Σ Ix², Sxy = Σ Ix·Iy, Syy = Σ Iy², Sxt = Σ Ix·It and Syt = Σ Iy·It, all summed over the window. A 2×2 system has a closed form, so there is no solver and no iteration — just a determinant:

det = Sxx·Syy − Sxy²
u   = (Sxy·Syt − Syy·Sxt) / det
v   = (Sxy·Sxt − Sxx·Syt) / det

Every pixel solves its own tiny system, reading only its own neighbourhood and writing only its own cell. That is the ideal GPU shape: 4,096 independent 2×2 solves with no coordination whatsoever — the gather formulation, exactly as in Convolution & Filters, with a little linear algebra at the end.

And a warning the input is built to deliver. det goes to zero where the window has nothing to say: a flat patch (no gradient at all) or a stretch of parallel edges (all gradients pointing the same way — task 2's problem, unchanged by a bigger window). Divide by it anyway and you get NaN or vectors hundreds of pixels long. Guard it.

25 equations walk into a 2x2 system; one of them comes out with an answer — A five by five window of pixels, each contributing one brightness-constancy equation, gathered into five running sums that fill a two by two matrix and a right-hand side, solved once for the centre pixel.

Goal: solve the 5×5 Lucas–Kanade system for every pixel — plane 0 is u, plane 1 is v — and return 0 wherever the determinant is too small to trust.

Requirements

Accumulate Sxx, Sxy, Syy, Sxt, Syt over the 5×5 window, sample indexes clamped to 0…this.constants.last
Compute det = Sxx * Syy - Sxy * Sxy
If Math.abs(det) < this.constants.eps, return 0 — the window has nothing to say
Otherwise plane 0 returns (Sxy * Syt - Syy * Sxt) / det and plane 1 returns (Sxy * Sxt - Sxx * Syt) / det
Same sign convention as task 2: (u, v) is where the content went, positive u rightward and positive v downward

Hint 1 — the window loop

The same clamped double loop the box blur used, five wide instead of three:

for (let wy = 0; wy < 5; wy++) {
  for (let wx = 0; wx < 5; wx++) {
    let sy = this.thread.y + wy - 2;
    if (sy < 0) sy = 0;
    if (sy > this.constants.last) sy = this.constants.last;
    // …same for sx, then read the three derivatives at [sy][sx]
  }
}

Hint 2 — five running sums, not two

Inside the loop, all five products accumulate together:

const ix = derivs[0][sy][sx];
const iy = derivs[1][sy][sx];
const it = derivs[2][sy][sx];
sxx += ix * ix;
sxy += ix * iy;
syy += iy * iy;
sxt += ix * it;
syt += iy * it;

sxy is the one that gets forgotten. It is the off-diagonal term — the thing that couples u and v — and dropping it silently turns the 2×2 solve into two unrelated divisions.

Hint 3 — the guard and the two returns

const det = sxx * syy - sxy * sxy;
if (Math.abs(det) < this.constants.eps) {
  return 0;
}
if (this.thread.z === 0) {
  return (sxy * syt - syy * sxt) / det;
}
return (sxy * sxt - sxx * syt) / det;

Compute the sums once, before the branch on z — both components need all five of them.

Same idea elsewhere

This is the shape GPUs were built for: a fixed-size neighbourhood read, a handful of registers, a closed-form solve, no synchronisation. A CUDA implementation stages the derivative tiles in shared memory and keeps the five sums in registers; a WGSL compute shader does the same with a workgroup var<workgroup> tile. gpu.js has neither, so every thread re-reads its own window — more traffic, identical answer, and the algorithm is unchanged.

All tasks in Optical Flow

This page is an interactive exercise — the editor, the GPU runner and your saved progress need JavaScript. The text above is the full brief.