RecipesCam/src/utils/burstAlign.ts

// Alignment maths for the low-light hand-held burst. Pure — no Skia, no React
// Native — so the check beside it runs on plain node:
//   node src/utils/burstAlign.check.ts
//
// Frames taken within a second of each other differ by translation only (the
// zoom and the lens do not change between them), so the estimate is one shift
// per axis, read off the row and column luminance PROFILES of a downscaled
// frame. Correlating two profiles costs O(w+h) per candidate offset against
// O(w*h) for 2D block matching — the difference between a few hundred ms and
// minutes of JS thread for a 12MP burst. The 2D check then confirms the profile
// estimate against the actual pixels, so a periodic scene (blinds, railings,
// a tiled wall) cannot lock the estimator onto the wrong repetition — and (0,0)
// stays a candidate throughout, so a scene neither stage can read merges
// unshifted rather than wrongly shifted.

export interface Profiles {
  /** Mean luminance of each row, top to bottom. */
  rows: Float32Array;
  /** Mean luminance of each column, left to right. */
  cols: Float32Array;
}

/** RGBA pixels as they come out of Skia's readPixels (8-bit or float). */
export type Pixels = Uint8Array | Float32Array;

// Rec.709 luma on the unmanaged 8-bit values. A monotone luminance proxy is all
// a correlation needs — this is not a colour-managed Y.
const LUMA_R = 0.2126;
const LUMA_G = 0.7152;
const LUMA_B = 0.0722;

export function profilesOf(rgba: Pixels, w: number, h: number): Profiles {
  const rows = new Float32Array(h);
  const cols = new Float32Array(w);
  for (let y = 0; y < h; y++) {
    const row = y * w * 4;
    let sum = 0;
    for (let x = 0; x < w; x++) {
      const p = row + x * 4;
      const l = rgba[p] * LUMA_R + rgba[p + 1] * LUMA_G + rgba[p + 2] * LUMA_B;
      sum += l;
      cols[x] += l;
    }
    rows[y] = sum / w;
  }
  for (let x = 0; x < w; x++) cols[x] /= h;
  return { rows, cols };
}

/**
 * The shift of `cur` against `ref`, in profile samples, searched over
 * `expect ± max` and scored by mean absolute difference over the overlap.
 * Mean, not sum, so a candidate that only matches on a handful of pixels is not
 * favoured; overlaps shorter than half the signal are rejected outright.
 * Convention: `ref[i] ≈ cur[i - shift]`.
 */
export function bestShift(ref: Float32Array, cur: Float32Array, max: number, expect = 0): number {
  const center = Math.round(expect);
  const lo = Math.max(-(cur.length - 1), center - max);
  const hi = Math.min(ref.length - 1, center + max);
  let best = center;
  let bestErr = Infinity;
  for (let s = lo; s <= hi; s++) {
    const from = Math.max(0, s);
    const to = Math.min(ref.length, cur.length + s);
    if (to - from < ref.length / 2) continue;
    let err = 0;
    for (let i = from; i < to; i++) err += Math.abs(ref[i] - cur[i - s]);
    err /= to - from;
    if (err < bestErr) {
      bestErr = err;
      best = s;
    }
  }
  return best;
}

/**
 * Mean absolute luminance difference between two downscaled frames at a
 * candidate shift, over the overlap — the same convention as `bestShift`
 * (`a[y][x] ≈ b[y - dy][x - dx]`). `Infinity` when the candidate keeps less than
 * half the frame: a nearly disjoint pair would otherwise score a suspiciously
 * low error on whatever flat area it did keep.
 */
export function sadAt(a: Pixels, b: Pixels, w: number, h: number, dx: number, dy: number): number {
  const x0 = Math.max(0, dx);
  const y0 = Math.max(0, dy);
  const x1 = Math.min(w, w + dx);
  const y1 = Math.min(h, h + dy);
  if (x1 <= x0 || y1 <= y0) return Infinity;
  if ((x1 - x0) * (y1 - y0) < (w * h) / 2) return Infinity;
  let sum = 0;
  for (let y = y0; y < y1; y++) {
    const ra = (y * w + x0) * 4;
    const rb = ((y - dy) * w + (x0 - dx)) * 4;
    for (let x = x0; x < x1; x++) {
      const pa = ra + (x - x0) * 4;
      const pb = rb + (x - x0) * 4;
      const la = a[pa] * LUMA_R + a[pa + 1] * LUMA_G + a[pa + 2] * LUMA_B;
      const lb = b[pb] * LUMA_R + b[pb + 1] * LUMA_G + b[pb + 2] * LUMA_B;
      sum += Math.abs(la - lb);
    }
  }
  return sum / ((x1 - x0) * (y1 - y0));
}

/**
 * Confirm a profile estimate on the pixels: search `±radius` around it, keeping
 * (0,0) in the running, and return the best 2D match. This is what makes a
 * profile lock onto a real translation rather than onto a repeated pattern —
 * and what makes "no shift" win whenever the pixels do not support a shift.
 */
export function nearest2D(
  ref: Pixels,
  cur: Pixels,
  w: number,
  h: number,
  seedX: number,
  seedY: number,
  radius: number
): { x: number; y: number } {
  let best = { x: 0, y: 0 };
  let bestErr = sadAt(ref, cur, w, h, 0, 0);
  for (let dy = seedY - radius; dy <= seedY + radius; dy++) {
    for (let dx = seedX - radius; dx <= seedX + radius; dx++) {
      if (dx === 0 && dy === 0) continue;
      const err = sadAt(ref, cur, w, h, dx, dy);
      if (err < bestErr) {
        bestErr = err;
        best = { x: dx, y: dy };
      }
    }
  }
  return best;
}