calibrate: find the chart beside a shadow, and hold the matrix to the camera's noise

dist/calibrate.js

@@ -468,13 +468,49 @@ function summarise(errs, k) {
 	};
 }
 
+/*
+ * How much a live matrix amplifies the sensor's noise, per output channel.
+ *
+ * Each output is a weighted sum of the three white-balanced inputs, so its
+ * noise grows with the length of that row of weights -- with the white
+ * balance gains folded in, because they multiply the noise before the matrix
+ * sees it. A row of [1, 0, 0] passes red's noise through unchanged; a row
+ * that subtracts its neighbours to buy saturation amplifies all three.
+ *
+ * This is the cost a chart fit never sees. The 24 patches are averages, so
+ * their noise is gone before the solver looks at them, and the fit is free to
+ * buy accuracy with amplification. On a lab gk7605v100 + SC2239 under a
+ * 2330 K lamp it did: the fit took mean ΔE2000 from 8.9 to 5.7 with a green
+ * row that amplified 3.73 against the 1.69 the camera ran before, and the
+ * camera's video carried 2.2 to 2.6 times the chroma noise on flat surfaces,
+ * measured over 40 frames against the matrix it replaced.
+ */
+export function noiseGain(ccm, neutral) {
+	const g = [1 / neutral[0], 1 / neutral[1], 1 / neutral[2]];
+	return [0, 1, 2].map((r) =>
+		Math.hypot(ccm[r * 3] * g[0], ccm[r * 3 + 1] * g[1], ccm[r * 3 + 2] * g[2]));
+}
+
+/*
+ * How much noisier than the camera's own matrix a calibration may make the
+ * picture, by default. Measured on the same camera, same light, same chart:
+ * held at 1.25 times, video chroma noise rose 0 to 30 per cent instead of 2.2
+ * to 2.6 times, and the chart's colour error came out at 10.4 against 10.8
+ * for the unheld matrix -- the amplification bought almost nothing the
+ * picture shows. Held at 1.0 the error was 12.9, and the camera's own 14.0.
+ */
+export const NOISE_HEADROOM = 1.25;
+
 /*
  * measured: 24 camera RGB triples, black-subtracted, in chart order.
  * opts.clipped: per patch, the fraction of it the sampler found clipped.
  * opts.colorMatrices: the DNG's ColorMatrix1/2 with their illuminants, so the
  *   light can be named.
  * opts.cct: the light's temperature, when someone knows it better than the
  *   camera does.
+ * opts.noiseBudget: the most each output channel may amplify the sensor's
+ *   noise, as noiseGain() measures it. The fit is held under it, and what the
+ *   chart alone would have asked for comes back as `noise.free`.
  *
  * Returns the white balance the neutral row implies, both matrices, the light,
  * and how well the result actually fits -- because a solve always returns
@@ -575,6 +611,54 @@ export function solveFromPatches(measured, opts = {}) {
 		throw new Error(`these patches fit the chart at a mean of ${fit.meanDeltaE.toFixed(1)} ΔE2000, ` +
 			`and a chart fits under ${MAX_MEAN_DELTA_E} — the corners are not on a colour chart`);
 
+	/*
+	 * Held under the noise budget, when there is one. The chart has been
+	 * judged above on the free fit -- whether these are a chart's patches does
+	 * not depend on how much noise the answer may cost.
+	 *
+	 * The budget is a penalty on each row's excess, made steeper until it
+	 * holds: a fixed weight left rows a few per cent over, and "a few per cent
+	 * over" is not a limit. Started from the free fit, which is the nearest
+	 * answer and the one a slack budget returns unchanged.
+	 */
+	let ccmOut = ccm, fitOut = fit, noise = null;
+	const free = noiseGain(ccm, neutral);
+	if (opts.noiseBudget) {
+		const budget = opts.noiseBudget;
+		if (budget.length !== 3 || budget.some((v) => !(v > 0)))
+			throw new Error('a noise budget is three positive numbers');
+		let q = p;
+		for (let mu = 30; free.some((v, r) => v > budget[r]) && mu <= 3e5; mu *= 10) {
+			const held = (x) => residual(x).concat(noiseGain(ccmFrom(x), neutral)
+				.map((v, r) => mu * Math.max(0, v - budget[r])));
+			q = levenberg(held, q);
+			if (noiseGain(ccmFrom(q), neutral).every((v, r) => v <= budget[r] * 1.005)) break;
+		}
+		/* Checked, not assumed: the steepest penalty can still leave a row
+		 * over a budget no row-sums-to-one matrix can meet, and a matrix
+		 * reported as held must be one. */
+		const reached = noiseGain(ccmFrom(q), neutral);
+		const over = Math.max(...reached.map((v, r) => v / budget[r]));
+		if (over > 1.005)
+			throw new Error(`no colour matrix fits this chart within the noise allowed — the ` +
+				`closest is ${over.toFixed(2)}× over it`);
+		if (q !== p) {
+			ccmOut = ccmFrom(q);
+			const errsHeld = balanced.map((_, i) => weights[i] ?
+				deltaE2000(lab(apply3(LINEAR_SRGB_TO_XYZ50, apply3(ccmOut, balanced[i])
+					.map((v) => v * Math.exp(q[6])))), CHART_LAB50[i]) : null);
+			fitOut = summarise(errsHeld, Math.exp(q[6]));
+			fitOut.patches = fit.patches;
+			fitOut.meanDeltaE76 = errsHeld.reduce((s, e, i) => s + (e === null ? 0 :
+				deltaE(apply3(LINEAR_SRGB_TO_XYZ50, apply3(ccmOut, balanced[i])
+					.map((v) => v * Math.exp(q[6]))), CHART_XYZ50[i])), 0) / fit.patches;
+		}
+		noise = { gain: noiseGain(ccmOut, neutral), budget: budget.slice(),
+			held: q !== p, free: q !== p ? { ccm, fit, gain: free } : null };
+	} else {
+		noise = { gain: free, budget: null, held: false, free: null };
+	}
+
 	/*
 	 * ColorMatrix1 maps XYZ *under the light the chart was lit by* to the
 	 * camera. The reference is D50, so it is carried to that light first --
@@ -611,9 +695,9 @@ export function solveFromPatches(measured, opts = {}) {
 	const colorMatrix = fitted.map((v) => v / peak);
 
 	return {
-		neutral, colorMatrix, ccm, balanced, weights,
+		neutral, colorMatrix, ccm: ccmOut, balanced, weights,
 		light, estimated,
-		fit,
+		fit: fitOut, noise,
 	};
 }
 

dist/editor.js

@@ -12,7 +12,7 @@
  * temperature/tint pair it has no model for yet.
  */
 
-import { solveFromPatches, patchCentres, scoreCcm, CHART_COLS, CHART_ROWS } from './calibrate.js';
+import { solveFromPatches, patchCentres, scoreCcm, noiseGain, NOISE_HEADROOM, CHART_COLS, CHART_ROWS } from './calibrate.js';
 import { parseIni, readColour, fitAwbCurve, gainsForCt, mergeCcmTables, colourFragment,
 	readDefectCorrection, enableDefectCorrection } from './iqprofile.js';
 import { summarise, peakHold, normalise, coarsen, sweepZones, zoneDetail,
@@ -2849,10 +2849,46 @@ export function mountEditor(root, {
 			patches.push(got.raw);
 			clipped.push(got.clipped || 0);
 		}
-		return solveFromPatches(patches, {
-			clipped,
-			colorMatrices: state.info && state.info.colorMatrices,
+		const opts = { clipped, colorMatrices: state.info && state.info.colorMatrices };
+		const free = solveFromPatches(patches, opts);
+		/*
+		 * No noisier than the camera is today, give or take NOISE_HEADROOM.
+		 *
+		 * The chart cannot say what a matrix costs in noise -- its patches are
+		 * averages -- so the fit is held against the matrix the camera already
+		 * runs at this light, read out of its own profile. A camera whose
+		 * profile cannot be read, or a light that cannot be named, is solved
+		 * free, and the result says so rather than implying it was checked.
+		 */
+		let own = null, why = null;
+		if (!(calibrate && calibrate.baseline)) why = 'no-baseline';
+		else if (!free.light) why = 'no-light';
+		else {
+			/* Only reading the camera's profile may fall back to an unheld fit.
+			 * A budget that cannot be met is an answer, and it is shown as one. */
+			try {
+				own = readColour(parseIni(await calibrate.baseline())).ccm;
+			} catch {
+				own = null;
+			}
+			if (!own) why = 'no-baseline';
+		}
+		if (!own) {
+			free.noise.unchecked = why;
+			return free;
+		}
+		/* The camera's matrix at this light, through the white balance the
+		 * camera itself chose for this frame -- its picture as it is, rather
+		 * than as it would be with the chart's balance. */
+		const at = vendorCcmAt(own, free.light.cct);
+		const shot = state.info && state.info.neutral;
+		const camNeutral = shot && shot.every((v) => v > 0) ? shot : free.neutral;
+		const camera = noiseGain(at, camNeutral);
+		const held = solveFromPatches(patches, {
+			...opts, noiseBudget: camera.map((v) => v * NOISE_HEADROOM),
 		});
+		held.noise.camera = camera;
+		return held;
 	}
 
 	/* The panel, and the way back.
@@ -3057,6 +3093,35 @@ export function mountEditor(root, {
 			out.append(lt);
 		}
 
+		/*
+		 * What the colour cost in noise, in words a camera's owner can act
+		 * on: how much noisier than today, and what matching the chart
+		 * exactly would have cost instead.
+		 */
+		{
+			const nz = solved.noise;
+			const worst = (g, ref) => Math.max(...g.map((v, k) => v / ref[k]));
+			const n = el('p', 're-note');
+			n.style.margin = '0 0 6px';
+			n.dataset.role = 'noise';
+			if (nz && nz.camera) {
+				const now = worst(nz.gain, nz.camera);
+				n.textContent = nz.held
+					? `Kept to ${now.toFixed(2)}× the noise of the camera's picture as it is now — ` +
+						'a calibration saved earlier counts as "now". ' +
+						`Matching the chart as closely as possible would fit at ΔE2000 ` +
+						`${nz.free.fit.meanDeltaE.toFixed(2)} and make the video ` +
+						`${worst(nz.free.gain, nz.camera).toFixed(1)}× as noisy.`
+					: `${now.toFixed(2)}× the noise of the camera's picture as it is now, ` +
+						`inside the ${NOISE_HEADROOM}× allowed.`;
+			} else {
+				n.textContent = 'Not checked against the camera\'s own colour, so how much ' +
+					'noisier this makes the picture is unknown: a close fit on a chart can ' +
+					'amplify the sensor\'s noise several times over.';
+			}
+			out.append(n);
+		}
+
 		out.append(Object.assign(el('h3', 're-cap'), { textContent: 'Live matrix, camera to display' }));
 		out.append(matrixTable(solved.ccm));
 		out.append(Object.assign(el('h3', 're-cap'), { textContent: 'ColorMatrix1, XYZ to camera' }));

src/calibrate.js

@@ -468,13 +468,49 @@ function summarise(errs, k) {
 	};
 }
 
+/*
+ * How much a live matrix amplifies the sensor's noise, per output channel.
+ *
+ * Each output is a weighted sum of the three white-balanced inputs, so its
+ * noise grows with the length of that row of weights -- with the white
+ * balance gains folded in, because they multiply the noise before the matrix
+ * sees it. A row of [1, 0, 0] passes red's noise through unchanged; a row
+ * that subtracts its neighbours to buy saturation amplifies all three.
+ *
+ * This is the cost a chart fit never sees. The 24 patches are averages, so
+ * their noise is gone before the solver looks at them, and the fit is free to
+ * buy accuracy with amplification. On a lab gk7605v100 + SC2239 under a
+ * 2330 K lamp it did: the fit took mean ΔE2000 from 8.9 to 5.7 with a green
+ * row that amplified 3.73 against the 1.69 the camera ran before, and the
+ * camera's video carried 2.2 to 2.6 times the chroma noise on flat surfaces,
+ * measured over 40 frames against the matrix it replaced.
+ */
+export function noiseGain(ccm, neutral) {
+	const g = [1 / neutral[0], 1 / neutral[1], 1 / neutral[2]];
+	return [0, 1, 2].map((r) =>
+		Math.hypot(ccm[r * 3] * g[0], ccm[r * 3 + 1] * g[1], ccm[r * 3 + 2] * g[2]));
+}
+
+/*
+ * How much noisier than the camera's own matrix a calibration may make the
+ * picture, by default. Measured on the same camera, same light, same chart:
+ * held at 1.25 times, video chroma noise rose 0 to 30 per cent instead of 2.2
+ * to 2.6 times, and the chart's colour error came out at 10.4 against 10.8
+ * for the unheld matrix -- the amplification bought almost nothing the
+ * picture shows. Held at 1.0 the error was 12.9, and the camera's own 14.0.
+ */
+export const NOISE_HEADROOM = 1.25;
+
 /*
  * measured: 24 camera RGB triples, black-subtracted, in chart order.
  * opts.clipped: per patch, the fraction of it the sampler found clipped.
  * opts.colorMatrices: the DNG's ColorMatrix1/2 with their illuminants, so the
  *   light can be named.
  * opts.cct: the light's temperature, when someone knows it better than the
  *   camera does.
+ * opts.noiseBudget: the most each output channel may amplify the sensor's
+ *   noise, as noiseGain() measures it. The fit is held under it, and what the
+ *   chart alone would have asked for comes back as `noise.free`.
  *
  * Returns the white balance the neutral row implies, both matrices, the light,
  * and how well the result actually fits -- because a solve always returns
@@ -575,6 +611,54 @@ export function solveFromPatches(measured, opts = {}) {
 		throw new Error(`these patches fit the chart at a mean of ${fit.meanDeltaE.toFixed(1)} ΔE2000, ` +
 			`and a chart fits under ${MAX_MEAN_DELTA_E} — the corners are not on a colour chart`);
 
+	/*
+	 * Held under the noise budget, when there is one. The chart has been
+	 * judged above on the free fit -- whether these are a chart's patches does
+	 * not depend on how much noise the answer may cost.
+	 *
+	 * The budget is a penalty on each row's excess, made steeper until it
+	 * holds: a fixed weight left rows a few per cent over, and "a few per cent
+	 * over" is not a limit. Started from the free fit, which is the nearest
+	 * answer and the one a slack budget returns unchanged.
+	 */
+	let ccmOut = ccm, fitOut = fit, noise = null;
+	const free = noiseGain(ccm, neutral);
+	if (opts.noiseBudget) {
+		const budget = opts.noiseBudget;
+		if (budget.length !== 3 || budget.some((v) => !(v > 0)))
+			throw new Error('a noise budget is three positive numbers');
+		let q = p;
+		for (let mu = 30; free.some((v, r) => v > budget[r]) && mu <= 3e5; mu *= 10) {
+			const held = (x) => residual(x).concat(noiseGain(ccmFrom(x), neutral)
+				.map((v, r) => mu * Math.max(0, v - budget[r])));
+			q = levenberg(held, q);
+			if (noiseGain(ccmFrom(q), neutral).every((v, r) => v <= budget[r] * 1.005)) break;
+		}
+		/* Checked, not assumed: the steepest penalty can still leave a row
+		 * over a budget no row-sums-to-one matrix can meet, and a matrix
+		 * reported as held must be one. */
+		const reached = noiseGain(ccmFrom(q), neutral);
+		const over = Math.max(...reached.map((v, r) => v / budget[r]));
+		if (over > 1.005)
+			throw new Error(`no colour matrix fits this chart within the noise allowed — the ` +
+				`closest is ${over.toFixed(2)}× over it`);
+		if (q !== p) {
+			ccmOut = ccmFrom(q);
+			const errsHeld = balanced.map((_, i) => weights[i] ?
+				deltaE2000(lab(apply3(LINEAR_SRGB_TO_XYZ50, apply3(ccmOut, balanced[i])
+					.map((v) => v * Math.exp(q[6])))), CHART_LAB50[i]) : null);
+			fitOut = summarise(errsHeld, Math.exp(q[6]));
+			fitOut.patches = fit.patches;
+			fitOut.meanDeltaE76 = errsHeld.reduce((s, e, i) => s + (e === null ? 0 :
+				deltaE(apply3(LINEAR_SRGB_TO_XYZ50, apply3(ccmOut, balanced[i])
+					.map((v) => v * Math.exp(q[6]))), CHART_XYZ50[i])), 0) / fit.patches;
+		}
+		noise = { gain: noiseGain(ccmOut, neutral), budget: budget.slice(),
+			held: q !== p, free: q !== p ? { ccm, fit, gain: free } : null };
+	} else {
+		noise = { gain: free, budget: null, held: false, free: null };
+	}
+
 	/*
 	 * ColorMatrix1 maps XYZ *under the light the chart was lit by* to the
 	 * camera. The reference is D50, so it is carried to that light first --
@@ -611,9 +695,9 @@ export function solveFromPatches(measured, opts = {}) {
 	const colorMatrix = fitted.map((v) => v / peak);
 
 	return {
-		neutral, colorMatrix, ccm, balanced, weights,
+		neutral, colorMatrix, ccm: ccmOut, balanced, weights,
 		light, estimated,
-		fit,
+		fit: fitOut, noise,
 	};
 }