|
| 1 | +// Pigment (subtractive) mixing via Kubelka–Munk theory on derived spectra. |
| 2 | +// |
| 3 | +// Two pieces, both derivation rather than tables: |
| 4 | +// |
| 5 | +// 1. Spectral upsampling (Jakob & Hanika 2019, "A Low-Dimensional Function |
| 6 | +// Space for Efficient Spectral Upsampling"): a color's plausible |
| 7 | +// reflectance is R(λ) = sigmoid(c₂ + c₁·λ̃ + c₀·λ̃²), with the three |
| 8 | +// coefficients NEWTON-SOLVED per color so that integrating R under D65 |
| 9 | +// against our CMFs reproduces the color's XYZ exactly. The paper |
| 10 | +// precomputes lookup tables; we solve directly (analytic Jacobian, |
| 11 | +// damped Newton) — no transcribed coefficients anywhere. |
| 12 | +// |
| 13 | +// 2. Kubelka–Munk single-constant mixing (Kubelka & Munk 1931; Duncan 1940 |
| 14 | +// for mixtures): per wavelength, K/S = (1−R)²/(2R); concentrations mix |
| 15 | +// linearly in K/S; back via R = 1 + K/S − √((K/S)² + 2·K/S). |
| 16 | +// |
| 17 | +// `t` is pigment concentration, not perceptual position — that's the |
| 18 | +// physics. Mixing yellow and blue passes through green, as paint does. |
| 19 | + |
| 20 | +import { CMF_1931_2, D65_SPD } from './data.js'; |
| 21 | +import { sampleSpd, reflectanceToXyz } from './index.js'; |
| 22 | +import { convert } from '../core/convert.js'; |
| 23 | +import { invert, mulVec } from '../core/mat3.js'; |
| 24 | + |
| 25 | +const N = CMF_1931_2.x.length; |
| 26 | +const START = CMF_1931_2.start; |
| 27 | +const STEP = CMF_1931_2.step; |
| 28 | + |
| 29 | +// Precomputed per-grid-point: normalized wavelength basis and D65-weighted |
| 30 | +// CMF products with the perfect-reflector normalization folded in. |
| 31 | +const LT = new Array(N); // λ̃ ∈ [-1, 1], centered for conditioning |
| 32 | +const SX = new Array(N); |
| 33 | +const SY = new Array(N); |
| 34 | +const SZ = new Array(N); |
| 35 | +{ |
| 36 | + let k = 0; |
| 37 | + for (let i = 0; i < N; i++) { |
| 38 | + const lambda = START + i * STEP; |
| 39 | + LT[i] = (lambda - 595) / 235; |
| 40 | + const s = sampleSpd(D65_SPD, lambda); |
| 41 | + SX[i] = s * CMF_1931_2.x[i]; |
| 42 | + SY[i] = s * CMF_1931_2.y[i]; |
| 43 | + SZ[i] = s * CMF_1931_2.z[i]; |
| 44 | + k += SY[i]; |
| 45 | + } |
| 46 | + for (let i = 0; i < N; i++) { |
| 47 | + SX[i] /= k; SY[i] /= k; SZ[i] /= k; |
| 48 | + } |
| 49 | +} |
| 50 | + |
| 51 | +const sigmoid = (x) => 0.5 + x / (2 * Math.sqrt(1 + x * x)); |
| 52 | +const sigmoidInv = (y) => { |
| 53 | + const u = 2 * y - 1; |
| 54 | + return u / Math.sqrt(1 - u * u); |
| 55 | +}; |
| 56 | + |
| 57 | +// Integrate the sigmoid-polynomial reflectance and its Jacobian. |
| 58 | +function evalF(c, F, J) { |
| 59 | + F[0] = 0; F[1] = 0; F[2] = 0; |
| 60 | + for (let r = 0; r < 9; r++) J[r] = 0; |
| 61 | + for (let i = 0; i < N; i++) { |
| 62 | + const lt = LT[i]; |
| 63 | + const p = c[0] * lt * lt + c[1] * lt + c[2]; |
| 64 | + const R = sigmoid(p); |
| 65 | + F[0] += R * SX[i]; F[1] += R * SY[i]; F[2] += R * SZ[i]; |
| 66 | + const d = 0.5 / Math.pow(1 + p * p, 1.5); // S'(p) |
| 67 | + const b0 = d * lt * lt, b1 = d * lt, b2 = d; |
| 68 | + J[0] += b0 * SX[i]; J[1] += b1 * SX[i]; J[2] += b2 * SX[i]; |
| 69 | + J[3] += b0 * SY[i]; J[4] += b1 * SY[i]; J[5] += b2 * SY[i]; |
| 70 | + J[6] += b0 * SZ[i]; J[7] += b1 * SZ[i]; J[8] += b2 * SZ[i]; |
| 71 | + } |
| 72 | +} |
| 73 | + |
| 74 | +const _F = [0, 0, 0]; |
| 75 | +const _J = new Array(9); |
| 76 | +const _step = [0, 0, 0]; |
| 77 | +const _resid = [0, 0, 0]; |
| 78 | + |
| 79 | +function solveCoeffs(xyz) { |
| 80 | + // clamp luminance into the solvable open interval, preserving chromaticity |
| 81 | + let [X, Y, Z] = xyz; |
| 82 | + if (Y <= 1e-6) return [0, 0, sigmoidInv(1e-4)]; // black: flat dark reflectance |
| 83 | + const yClamp = Math.min(Math.max(Y, 5e-4), 0.9995); |
| 84 | + const s = yClamp / Y; |
| 85 | + X *= s; Y *= s; Z *= s; |
| 86 | + |
| 87 | + const c = [0, 0, sigmoidInv(Math.min(Math.max(Y, 1e-3), 0.999))]; |
| 88 | + let err = Infinity; |
| 89 | + for (let iter = 0; iter < 60; iter++) { |
| 90 | + evalF(c, _F, _J); |
| 91 | + _resid[0] = _F[0] - X; _resid[1] = _F[1] - Y; _resid[2] = _F[2] - Z; |
| 92 | + err = Math.max(Math.abs(_resid[0]), Math.abs(_resid[1]), Math.abs(_resid[2])); |
| 93 | + if (err < 1e-10) break; |
| 94 | + mulVec(invert(_J), _resid, _step); |
| 95 | + // damped: halve the step while it increases the residual |
| 96 | + let scale = 1; |
| 97 | + for (let h = 0; h < 10; h++) { |
| 98 | + const trial = [c[0] - scale * _step[0], c[1] - scale * _step[1], c[2] - scale * _step[2]]; |
| 99 | + evalF(trial, _F, _J); |
| 100 | + const e2 = Math.max(Math.abs(_F[0] - X), Math.abs(_F[1] - Y), Math.abs(_F[2] - Z)); |
| 101 | + if (e2 < err) { |
| 102 | + c[0] = trial[0]; c[1] = trial[1]; c[2] = trial[2]; |
| 103 | + break; |
| 104 | + } |
| 105 | + scale *= 0.5; |
| 106 | + if (h === 9) { c[0] -= scale * _step[0]; c[1] -= scale * _step[1]; c[2] -= scale * _step[2]; } |
| 107 | + } |
| 108 | + } |
| 109 | + return c; |
| 110 | +} |
| 111 | + |
| 112 | +const _cache = new Map(); |
| 113 | + |
| 114 | +/** |
| 115 | + * A plausible smooth reflectance spectrum for a color (Jakob–Hanika sigmoid, |
| 116 | + * Newton-solved so it integrates back to the color's XYZ under D65 exactly). |
| 117 | + * Cached per color. |
| 118 | + */ |
| 119 | +export function reflectanceOf(coords, space = 'srgb') { |
| 120 | + const xyz = convert(coords, space, 'xyz-d65'); |
| 121 | + const key = `${xyz[0].toFixed(6)},${xyz[1].toFixed(6)},${xyz[2].toFixed(6)}`; |
| 122 | + let spd = _cache.get(key); |
| 123 | + if (spd) return spd; |
| 124 | + const c = solveCoeffs(xyz); |
| 125 | + const values = new Array(N); |
| 126 | + for (let i = 0; i < N; i++) { |
| 127 | + const lt = LT[i]; |
| 128 | + values[i] = sigmoid(c[0] * lt * lt + c[1] * lt + c[2]); |
| 129 | + } |
| 130 | + spd = { start: START, step: STEP, values }; |
| 131 | + if (_cache.size > 512) _cache.clear(); |
| 132 | + _cache.set(key, spd); |
| 133 | + return spd; |
| 134 | +} |
| 135 | + |
| 136 | +/** |
| 137 | + * Kubelka–Munk mix of two reflectance spectra at concentration t |
| 138 | + * (0 → all a, 1 → all b). Spectra must share a grid (reflectanceOf's do). |
| 139 | + */ |
| 140 | +export function kmMixReflectance(ra, rb, t, out) { |
| 141 | + const n = ra.values.length; |
| 142 | + const values = out?.values ?? new Array(n); |
| 143 | + for (let i = 0; i < n; i++) { |
| 144 | + const a = Math.min(Math.max(ra.values[i], 1e-5), 1 - 1e-5); |
| 145 | + const b = Math.min(Math.max(rb.values[i], 1e-5), 1 - 1e-5); |
| 146 | + const ksA = ((1 - a) * (1 - a)) / (2 * a); |
| 147 | + const ksB = ((1 - b) * (1 - b)) / (2 * b); |
| 148 | + const ks = ksA + t * (ksB - ksA); |
| 149 | + values[i] = 1 + ks - Math.sqrt(ks * ks + 2 * ks); |
| 150 | + } |
| 151 | + return { start: ra.start, step: ra.step, values }; |
| 152 | +} |
| 153 | + |
| 154 | +const _xyzOut = [0, 0, 0]; |
| 155 | + |
| 156 | +/** |
| 157 | + * Mix two colors as PIGMENTS (subtractive, Kubelka–Munk) rather than as |
| 158 | + * light. Yellow + blue = green. Input/output in `space`; t is pigment |
| 159 | + * concentration of b. |
| 160 | + */ |
| 161 | +export function pigmentMix(a, b, t, space = 'srgb', out = [0, 0, 0]) { |
| 162 | + const ra = reflectanceOf(a, space); |
| 163 | + const rb = reflectanceOf(b, space); |
| 164 | + const rm = kmMixReflectance(ra, rb, t); |
| 165 | + reflectanceToXyz(rm, undefined, _xyzOut); |
| 166 | + return convert(_xyzOut, 'xyz-d65', space, out); |
| 167 | +} |
0 commit comments