# Lucas–Kanade: Buy a Second Equation

*Task 3 of 5 · [Optical Flow](https://gpu.rocks/learn/optical-flow-e85c6dfa.md) · GPU.js Learn*

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:

```js
| 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:

```js
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.

## Figures

- **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

**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:

```js
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:

```js
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

```js
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.

## Starter code

```js
// One 2×2 least-squares solve per pixel, from a 5×5 window of equations.
const gpu = new GPU({ mode });

const lucasKanade = gpu.createKernel(function (derivs) {
  const x = this.thread.x;
  const y = this.thread.y;

  let sxx = 0;
  let sxy = 0;
  let syy = 0;
  let sxt = 0;
  let syt = 0;

  // TODO: loop over the 5×5 window around (x, y), clamping both sample
  // coordinates to 0…this.constants.last, and accumulate the five sums.

  // TODO: det = sxx * syy - sxy * sxy;
  //       return 0 when Math.abs(det) < this.constants.eps,
  //       otherwise plane 0 → (sxy * syt - syy * sxt) / det
  //                 plane 1 → (sxy * sxt - sxx * syt) / det
  return 0;
}, {
  output: [64, 64, 2],
  constants: { last: 63, eps: 1e-6 },
});

const flow = await lucasKanade(derivs);
console.log('textured region, true motion is (1, 1):', flow[0][30][45], flow[1][30][45]);
console.log('flat region:', flow[0][10][6], flow[1][10][6]);
```

---

Interactive version: https://gpu.rocks/learn/optical-flow-e85c6dfa/3

[Previous task](https://gpu.rocks/learn/optical-flow-e85c6dfa/2.md) · [Next task](https://gpu.rocks/learn/optical-flow-e85c6dfa/4.md)