// The colour the camera itself wrote into the file.
//
// LibRaw hands back the sensor's own data and nothing about how the body chose to
// render it: the white balance it locked, its picture style, its tone curve. The
// file carries the answer anyway — the preview saved inside it, which is the
// frame the photographer saw and the one a desktop viewer shows. So the develop
// is matched to that preview: develop a grid of the frame, take the same grid out
// of the embedded preview, and fit the curve between them — one cubic a channel,
// on the sRGB-encoded value the develop leaves.
//
// A 3x3 was fitted ahead of that curve for a long time, least squared from one
// grid to the other, and it is what put the cast the studio was reported with:
// a develop carries its luminance and its colour in the same three numbers, the
// luminance is two orders of magnitude the larger of the two, and a least squares
// that is dominated by it pays for the luminance with the colour — it holds the
// channel means and the bright/dark axis while it collapses the axis the eye
// reads as colour. The camera's own picture style is a curve, not a matrix.
//
// Measured on the Olympus ORF this was reported on, standard deviation of R-G
// over the frame (the file's own preview 8.3, the camera's own JPEG 12.3, and an
// olive frame is one where this has gone to nothing):
//   develop as the sensor left it                       11.6
//   the curve alone, no 3x3                             11.9
//   the curve under the fitted 3x3 (what used to ship)   4.6
// and on the frame's own means, curve alone against the preview: luminance 156.6
// against 156.8, R/G 1.023 against 1.014, B/G 0.853 against 0.843.
//
// The fit is on the sRGB-encoded values, exposure included: the reference is the
// file's own frame, so matching its brightness is part of matching its colour.
// Fitted in that encoded space, not in linear light, because a least squares in
// linear light weighs the highlights and comes back with the midtones wrong —
// the same frames scored dE00 6.4 fitted encoded against 9.2 fitted linear. Any
// body, any picture style, and no fitted table to go stale — a preview is the one
// profile every RAW already carries.

// The grid the fit reads, and the size of the preview map handed to it.
export const MATCH_GRID = 64;

// A block at or past this has nothing left to say about a rendering.
const CLIP = 0.97;
const FLOOR = 0.03;

// What the shader draws when the frame has no curve the fit can trust.
export const FLAT_TONE = new Float32Array([0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0]);

// The curve a channel needs to reach the camera's, least squared: the develop's
// own value for that channel against the block the camera wrote there.
//
// The cubic is pinned at BOTH black (x=0 -> y=0) and white (x=1 -> y=1).
// Pinned at both ends, a blown pixel (1,1,1) stays pure white, and a deep shadow
// pixel (0,0,0) stays pure black without shifting to purple or magenta noise.
function basis(i: number, x: number): number {
  const x3 = x * x * x;
  return i === 0 ? x - x3 : x * x - x3;
}

function solve2x2(A: number[][], d: number[]): number[] {
  const det = A[0][0] * A[1][1] - A[0][1] * A[1][0] + 1e-9;
  const a = (d[0] * A[1][1] - d[1] * A[0][1]) / det;
  const b = (A[0][0] * d[1] - A[1][0] * d[0]) / det;
  return [a, b];
}

// The curve itself, in the shader's own order: r a b c d, then g, then b.
function fitTone(blocks: { x: number[]; t: number[] }[]): Float32Array {
  const out = new Float32Array(12);
  for (let c = 0; c < 3; c++) {
    const A = [[0, 0], [0, 0]];
    const d = [0, 0];
    for (const { x, t } of blocks) {
      const v = Math.min(1, Math.max(0, x[c]));
      const b = [basis(0, v), basis(1, v)];
      for (let p = 0; p < 2; p++) {
        for (let k = 0; k < 2; k++) A[p][k] += b[p] * b[k];
        d[p] += b[p] * (t[c] - v * v * v);
      }
    }
    const [a, b] = solve2x2(A, d);
    out.set([0, a, b, 1 - a - b], c * 4);
  }
  return out;
}

// A curve that comes back folded, or wandered far from the value it was handed,
// is a fit off blocks that fought each other, not a rendering: the caller opens
// the frame as the sensor left it instead. The floor sits below the zero a curve
// is allowed to reach — the develop clamps there anyway — and the slope is only
// asked not to turn back.
function curveIsSane(tone: Float32Array): boolean {
  for (const v of tone) if (!Number.isFinite(v) || Math.abs(v) > 4) return false;
  for (let c = 0; c < 3; c++) {
    const [a, b, cc, dd] = [tone[c * 4], tone[c * 4 + 1], tone[c * 4 + 2], tone[c * 4 + 3]];
    let prev = a;
    for (let i = 1; i <= 16; i++) {
      const x = i / 16;
      const y = a + x * (b + x * (cc + x * dd));
      if (y < prev - 0.12 || y < -0.15 || y > 1.15) return false;
      prev = y;
    }
  }
  return true;
}

// `dev` and `ref` are n x n RGBA maps of the same frame, one byte per channel,
// both sRGB-encoded by the same resampler. Null when the two cannot be related —
// too few blocks carry colour — and the caller develops the frame as the sensor
// left it.
export function toneMatch(dev: Uint8Array, ref: Uint8Array, n = MATCH_GRID): Float32Array | null {
  const blocks: { x: number[]; t: number[] }[] = [];
  let used = 0;
  for (let i = 0; i < n * n; i++) {
    const x = [dev[i * 4] / 255, dev[i * 4 + 1] / 255, dev[i * 4 + 2] / 255];
    const t = [ref[i * 4] / 255, ref[i * 4 + 1] / 255, ref[i * 4 + 2] / 255];
    // A block the camera clipped has no colour left to copy, one at the floor is
    // only the pedestal, and one the develop clipped has lost the ratio between
    // its channels — none of them say anything about the rendering.
    if (t[0] > CLIP && t[1] > CLIP && t[2] > CLIP) continue;
    if (x[0] > CLIP && x[1] > CLIP && x[2] > CLIP) continue;
    if (Math.max(t[0], t[1], t[2]) < FLOOR) continue;
    used++;
    blocks.push({ x, t });
  }
  // A fit off a handful of blocks is a fit off noise.
  if (used < (n * n) / 8) return null;
  const tone = fitTone(blocks);
  return curveIsSane(tone) ? tone : null;
}