Metamerism: two fabrics that match in the shop and not in daylight (Python)

Play: switch the lamp over two gray swatches

A gray jacket and gray trousers, bought together, match exactly under the shop's warm lamp, an incandescent bulb at 2,856 K, and outside in daylight the trousers turn greenish. Both swatches in the game are measured fabrics, the trousers changed just enough to match the jacket under the bulb. Tap a lamp, warm bulb, white LED, or daylight (keys 1, 2, 3), and choose who looks, the eye or a camera (key V). Under each swatch is the light it reflects, drawn thick, with the other swatch's thin and gray. The number is ΔE00, the color difference between the swatches: 0.00 under the bulb and 4.98 in daylight, with a verdict beside it, same below 1, barely up to 2, noticeable above. Yet the two curves differ under every lamp, so how can they give one color?

Download the game

Show code

// game.js: two gray swatches, a jacket and a pair of trousers, under three lamps, seen by the eye or
// by a camera. The notebook writes the spectra and the 3 x 3 matrices into the #data element; every
// color and number on screen is computed from them here, whenever a button changes the state.
"use strict";

const DATA = JSON.parse(document.getElementById("data").textContent);
const LAMPS = ["A", "LED-B3", "D65"];               // the buttons and the keys 1, 2, 3, in this order
const SWATCHES = ["jacket", "trousers"];
const RAD = Math.PI / 180;

// ---- the state: which lamp is on and who looks
const DEFAULTS = {lamp: "A", viewer: "eye"};
const state = {...DEFAULTS};

// ---- from a spectrum to three numbers: X = k Σ S R x̄ Δλ, and likewise Y and Z
function tristimulus(R, S, sens) {
  let norm = 0;
  for (let i = 0; i < S.length; i++) norm += S[i] * DATA.cmf[i][1];
  const k = 100 / norm;                             // a perfect white reflector gets Y = 100
  const out = [0, 0, 0];
  for (let i = 0; i < S.length; i++) {
    for (let j = 0; j < 3; j++) out[j] += k * S[i] * R[i] * sens[i][j];
  }
  return out;
}

const mul = (M, v) => M.map((row) => row[0] * v[0] + row[1] * v[1] + row[2] * v[2]);

// The eye's XYZ directly; the camera's R, G, B through its 3 x 3 matrix fitted under this lamp.
function seenXYZ(R, s) {
  const S = DATA.lamps[s.lamp];
  if (s.viewer === "eye") return tristimulus(R, S, DATA.cmf);
  return mul(DATA.camera_to_xyz[s.lamp], tristimulus(R, S, DATA.camera));
}

function lab(XYZ, white) {                          // CIELAB against the lamp's own white
  const f = (t) => (t > 216 / 24389 ? Math.cbrt(t) : (24389 / 27 * t + 16) / 116);
  const [fx, fy, fz] = XYZ.map((v, j) => f(v / white[j]));
  return [116 * fy - 16, 500 * (fx - fy), 200 * (fy - fz)];
}

// ---- CIEDE2000 (CIE 142-2001; implementation notes: Sharma, Wu, and Dalal 2005), angles in degrees
function deltaE00([L1, a1, b1], [L2, a2, b2]) {
  const cos = (deg) => Math.cos(deg * RAD), sin = (deg) => Math.sin(deg * RAD);
  const hue = (b, a) => (a === 0 && b === 0 ? 0 : (Math.atan2(b, a) / RAD + 360) % 360);
  const C7 = ((Math.hypot(a1, b1) + Math.hypot(a2, b2)) / 2) ** 7;
  const G = 0.5 * (1 - Math.sqrt(C7 / (C7 + 25 ** 7)));    // stretches a near the gray axis
  const C1 = Math.hypot((1 + G) * a1, b1), C2 = Math.hypot((1 + G) * a2, b2);
  const h1 = hue(b1, (1 + G) * a1), h2 = hue(b2, (1 + G) * a2);
  let dh = h2 - h1;                                 // the shorter way around the hue circle
  if (C1 * C2 === 0) dh = 0;
  else if (dh > 180) dh -= 360;
  else if (dh < -180) dh += 360;
  const dL = L2 - L1, dC = C2 - C1, dH = 2 * Math.sqrt(C1 * C2) * sin(dh / 2);
  const L = (L1 + L2) / 2, C = (C1 + C2) / 2;
  let h = h1 + h2;                                  // the mean hue, also the short way around
  if (C1 * C2 !== 0) h = (Math.abs(h1 - h2) > 180 ? h + (h < 360 ? 360 : -360) : h) / 2;
  const T = 1 - 0.17 * cos(h - 30) + 0.24 * cos(2 * h) + 0.32 * cos(3 * h + 6) - 0.2 * cos(4 * h - 63);
  const SL = 1 + (0.015 * (L - 50) ** 2) / Math.sqrt(20 + (L - 50) ** 2);
  const SC = 1 + 0.045 * C, SH = 1 + 0.015 * C * T;
  const RT = -2 * Math.sqrt(C ** 7 / (C ** 7 + 25 ** 7)) * sin(60 * Math.exp(-(((h - 275) / 25) ** 2)));
  return Math.sqrt((dL / SL) ** 2 + (dC / SC) ** 2 + (dH / SH) ** 2 + RT * (dC / SC) * (dH / SH));
}

// ---- to the screen: adapted from the lamp's white to daylight (Bradford), then sRGB, 8 bit
function toScreen(XYZ, lamp) {
  const linear = mul(DATA.to_srgb[lamp], XYZ.map((v) => v / 100));
  const clipped = linear.some((v) => v < 0 || v > 1);
  const rgb = linear.map((v) => {
    const c = Math.min(1, Math.max(0, v));
    return Math.round(255 * (c <= 0.0031308 ? 12.92 * c : 1.055 * c ** (1 / 2.4) - 0.055));
  });
  return {rgb, clipped};
}

function look(s) {                                  // everything the screen shows, from the state
  const S = DATA.lamps[s.lamp];
  const white = tristimulus(S.map(() => 1), S, DATA.cmf);
  const swatches = SWATCHES.map((name) => {
    const R = DATA.reflectance[name], XYZ = seenXYZ(R, s);
    return {name, light: S.map((v, i) => v * R[i]), lab: lab(XYZ, white), ...toScreen(XYZ, s.lamp)};
  });
  return {swatches, dE: deltaE00(swatches[0].lab, swatches[1].lab)};
}

const verdict = (dE) => (dE < 1 ? "same" : dE < 2 ? "barely" : "noticeable");

// ---- drawing
const canvas = document.getElementById("scene");
const ctx = canvas.getContext("2d");
const $ = (id) => document.getElementById(id);
let view = {};

function layout() {
  const W = document.documentElement.clientWidth, H = document.documentElement.clientHeight;
  const dpr = window.devicePixelRatio || 1;
  const h = Math.max(120, H - $("panel").offsetHeight);
  canvas.style.height = `${h}px`;
  canvas.width = Math.round(W * dpr);
  canvas.height = Math.round(h * dpr);
  view = {W, H: h, dpr, font: parseFloat(getComputedStyle(document.body).fontSize)};
  render();
}

function render() {
  const css = getComputedStyle(document.documentElement);
  const color = (name) => css.getPropertyValue(name).trim();
  const {W, H, font} = view, seen = look(state);
  ctx.setTransform(view.dpr, 0, 0, view.dpr, 0, 0);

  // the wall: the page's background, tinted faintly by the lamp as a camera without white balance sees it
  ctx.fillStyle = color("--bg");
  ctx.fillRect(0, 0, W, H);
  ctx.globalAlpha = 0.18;
  ctx.fillStyle = DATA.tint[state.lamp];
  ctx.fillRect(0, 0, W, H);
  ctx.globalAlpha = 1;

  const pad = Math.max(10, 0.03 * W), gap = pad;
  const colW = Math.min((W - 2 * pad - gap) / 2, 1.2 * H);
  const x0 = (W - 2 * colW - gap) / 2;
  const name = 1.5 * font, caption = 1.5 * font, ticks = 1.4 * font;   // text rows above and below
  const room = H - 2 * pad - name - caption - ticks;
  const swatchH = 0.62 * room, plotH = 0.38 * room;
  const top = pad + name, plotTop = top + swatchH + caption;
  const peak = Math.max(...seen.swatches.flatMap((w) => w.light));   // both curves on one scale
  const lam = DATA.lambda, n = lam.length;
  const px = (x, i) => x + (i / (n - 1)) * colW;

  ctx.font = `${font}px system-ui, sans-serif`;
  seen.swatches.forEach((w, k) => {
    const x = x0 + k * (colW + gap), stroke = color(k === 0 ? "--ink" : "--accent");
    ctx.fillStyle = color("--fg");
    ctx.textBaseline = "alphabetic";
    ctx.fillText(w.name, x, top - 0.45 * font);
    ctx.fillStyle = `rgb(${w.rgb.join(",")})`;      // the swatch: its computed screen color
    ctx.fillRect(x, top, colW, swatchH);
    ctx.strokeStyle = stroke;
    ctx.lineWidth = 2;
    ctx.strokeRect(x + 1, top + 1, colW - 2, swatchH - 2);

    // the light it reflects, S(λ) R(λ); the other swatch's light thin and gray, for comparison
    ctx.fillStyle = color("--muted");
    ctx.fillText("reflected light", x, plotTop - 0.45 * font);
    const y = (v) => plotTop + plotH * (1 - v / peak);
    const other = seen.swatches[1 - k].light;
    ctx.strokeStyle = color("--muted");
    ctx.beginPath();
    ctx.moveTo(x, plotTop + plotH);
    ctx.lineTo(x + colW, plotTop + plotH);
    ctx.stroke();
    for (const [curve, c, lw] of [[other, color("--muted"), 1], [w.light, stroke, 2.2]]) {
      ctx.strokeStyle = c;
      ctx.lineWidth = lw;
      ctx.beginPath();
      curve.forEach((v, i) => (i ? ctx.lineTo(px(x, i), y(v)) : ctx.moveTo(px(x, i), y(v))));
      ctx.stroke();
    }
    ctx.textBaseline = "top";
    ctx.textAlign = "center";
    for (const nm of [400, 700]) ctx.fillText(nm === 700 ? "700 nm" : "400", px(x, lam.indexOf(nm)), plotTop + plotH + 3);
    ctx.textAlign = "left";
  });

  // the readout and the pressed buttons
  $("de").textContent = seen.dE.toFixed(2);
  $("verdict").textContent = verdict(seen.dE);
  $("gamut").textContent = seen.swatches.some((w) => w.clipped) ? "out of screen gamut" : "";
  document.querySelectorAll("button[data-lamp]").forEach((b) => b.setAttribute("aria-pressed", b.dataset.lamp === state.lamp));
  document.querySelectorAll("button[data-viewer]").forEach((b) => b.setAttribute("aria-pressed", b.dataset.viewer === state.viewer));
}

// ---- controls and keys; nothing moves on its own, so the game redraws on a change only
function set(partial) {
  Object.assign(state, partial);
  render();
}
document.querySelectorAll("button[data-lamp]").forEach((b) => b.addEventListener("click", () => set({lamp: b.dataset.lamp})));
document.querySelectorAll("button[data-viewer]").forEach((b) => b.addEventListener("click", () => set({viewer: b.dataset.viewer})));
document.addEventListener("keydown", (e) => {
  const key = e.key.toLowerCase();
  if (["1", "2", "3"].includes(key)) set({lamp: LAMPS[Number(key) - 1]});
  else if (key === "v") set({viewer: state.viewer === "eye" ? "camera" : "eye"});
});

window.addEventListener("resize", layout);
window.matchMedia("(prefers-color-scheme: dark)").addEventListener("change", render);
layout();
window.game = {state, set, render};

A spectrum becomes three numbers

Switch the lamp and the curves under the swatches change shape: the warm bulb tilts the light toward the red, daylight spreads it evenly. Each swatch still shows one color. The stills use the game's own data, downloaded by a hidden cell from the CIE (International Commission on Illumination): the matching functions, the lamps, and 99 measured reflectances, 107 kB.

Show code

import json
import warnings
from pathlib import Path

import numpy as np
import matplotlib.pyplot as plt
warnings.filterwarnings("ignore", message="IProgress not found")    # tqdm, imported by pooch, wants ipywidgets
import pooch
import colour

plt.rcParams.update({
    "figure.figsize": (7, 3.6), "figure.dpi": 110,
    "axes.spines.top": False, "axes.spines.right": False,
    "axes.grid": True, "grid.alpha": 0.25,
    "font.size": 11, "lines.linewidth": 1.8,
})
INK, ACCENT, SECOND, MUTED = "#1f2a44", "#c8553d", "#2a7f9e", "#8a8f98"

# CIE data sets (CC BY-SA 4.0), pinned by the SHA-256 the CIE publishes with each file
CIE = "https://files.cie.co.at/Publications-datasets/"
FILES = {
    "cmf": ("CIE_xyz_1931_2deg.csv",                 # CIE 018:2019, DOI 10.25039/CIE.DS.xvudnb9b
            "fa663e3535a7e0763a745993a1f0a192eb0275ac46ad2d1befd7626841e713c1"),
    "A": ("CIE_std_illum_A_1nm.csv",                 # DOI 10.25039/CIE.DS.8jsxjrsn
          "61ef23fe146b8b665c74706717ab28cec7db6c9022993490bdc71991f43cb59b"),
    "D65": ("CIE_std_illum_D65.csv",                 # DOI 10.25039/CIE.DS.hjfjmt59
            "e76f210bffff3d552ef7113025da5f325d5dfec200dd4b878b1a2f3a507032cb"),
    "LED": ("CIE_illum_LEDs.csv",                    # CIE 015:2018, DOI 10.25039/CIE.DS.vgssnyfg
            "6147313bd1ec570ef308df222e1fe1e58115fde1ad0a7d4ba91fb9dcb776b267"),
    "samples": ("CIE_srf_cfi.csv",                   # CIE 224:2017, DOI 10.25039/CIE.DS.wi5idbqu
                "eb1d04a0d489918970f32354c2261581e9f9ba59456e0feae7227c7b5413aa3b"),
}
pooch.get_logger().setLevel("WARNING")
paths = {key: pooch.retrieve(CIE + name, known_hash=f"sha256:{sha}") for key, (name, sha) in FILES.items()}
table = {key: np.loadtxt(path, delimiter=",") for key, path in paths.items()}

lam = np.arange(380, 781, 5)                         # nm; the 1 nm tables are read at every fifth value

def on_grid(t, column=1):
    rows = np.searchsorted(t[:, 0], lam)
    assert np.array_equal(t[rows, 0], lam)
    return t[rows, column]

cmf = np.column_stack([on_grid(table["cmf"], j) for j in (1, 2, 3)])   # x̄, ȳ, z̄
lamps = {"A": on_grid(table["A"]), "LED-B3": on_grid(table["LED"], 3), "D65": on_grid(table["D65"])}
samples = table["samples"][:, 1:]                    # 99 reflectances, one per column
jacket, green = samples[:, 61], samples[:, 46]       # samples 62 and 47, both textiles

# The Nikon D5100 measured at NPL (Darrodi et al. 2015), as colour bundles it
nikon = colour.MSDS_CAMERA_SENSITIVITIES["Nikon 5100 (NPL)"]
assert np.array_equal(nikon.wavelengths, lam)
camera = nikon.values                                # 81 x 3: R, G, B

size = sum(Path(p).stat().st_size for p in paths.values())
print(f"{len(FILES)} CIE files, {size / 1000:.0f} kB; {lam.size} wavelengths, {samples.shape[1]} samples, "
      f"lamps {', '.join(lamps)}; camera curves {camera.shape}")

5 CIE files, 107 kB; 81 wavelengths, 99 samples, lamps A, LED-B3, D65; camera curves (81, 3)

Show code

S = lamps["A"]
k = 100 / (S * cmf[:, 1]).sum()
light = S * jacket                                   # the light reaching the eye
scale = 0.75 * cmf[:, 0].max() / (light * cmf[:, 0]).max()   # one factor for all three fills: areas stay comparable
fig, (top, bottom) = plt.subplots(2, 1, figsize=(7, 4.4), sharex=True)
top.plot(lam, jacket, color=INK)
top.plot(lam, S / S.max(), color=MUTED, lw=1.2)
top.text(595, 0.17, "jacket", color=INK, ha="center")
top.text(600, 0.78, "warm bulb, scaled", color=MUTED, ha="right")
top.set(ylabel="reflectance", ylim=(0, 1.05))
spots = [(628, 1.02, "left"), (560, 1.2, "center"), (478, 1.6, "left")]   # label positions, x̄ ȳ z̄
for j, (name, ls) in enumerate(zip(["x̄", "ȳ", "z̄"], ["-", "--", ":"])):
    bottom.plot(lam, cmf[:, j], color=INK, ls=ls, lw=1.2)
    bottom.fill_between(lam, scale * light * cmf[:, j], color=INK, alpha=0.18, lw=0)
    bottom.plot(lam, scale * light * cmf[:, j], color=INK, ls=ls, lw=0.8)   # its edge, styled like its function
    value = k * (light * cmf[:, j]).sum()
    x, y, ha = spots[j]
    bottom.text(x, y, f"{name}: {'XYZ'[j]} = {value:.2f}", color=INK, ha=ha)
bottom.set(xlabel="λ / nm", ylabel="matching function", xlim=(380, 780), ylim=(0, 2.05))
plt.show()

Top: the jacket's reflectance from 380 to 780 nm, 0.3 to 0.47 across the middle and rising beyond 680 nm, and the warm bulb's spectrum rising toward the red. Bottom: the matching functions x̄, ȳ, z̄, each with the jacket's light weighted by it shaded underneath. The three areas give X = 35.16, Y = 32.83, Z = 13.01.

The eye has three kinds of cones, so a color is three numbers. The CIE 1931 matching functions x̄, ȳ, z̄ are a standard observer's three responses, to a good approximation a linear recombination of the cone sensitivities. Each number weighs the light reaching the eye, lamp spectrum S(λ) times reflectance R(λ), with one matching function and sums over wavelength, the shaded areas in the still:

\[X = k\sum_\lambda S(\lambda)\,R(\lambda)\,\bar x(\lambda),\qquad k = \frac{100}{\sum_\lambda S(\lambda)\,\bar y(\lambda)},\]

and the same with ȳ for Y and z̄ for Z. The constant k scales each lamp so that a perfect white reflector has Y = 100. The three sums are one matrix acting on the reflectance, XYZ = M_S R, with 3 rows and 81 columns, one per 5 nm from 380 to 780 nm. In the code, lamp[:, None] multiplies each row of the matching functions by its lamp value (broadcasting, see Vectorizing loops with NumPy).

def tristimulus_matrix(lamp):
    """The 3 x 81 matrix that turns a reflectance into X, Y, Z under this lamp."""
    k = 100 / (lamp * cmf[:, 1]).sum()               # a perfect white reflector gets Y = 100
    return k * (cmf * lamp[:, None]).T

M_A = tristimulus_matrix(lamps["A"])
X, Y, Z = M_A @ jacket
print(f"M_A: {M_A.shape[0]} x {M_A.shape[1]}")
print(f"jacket under the bulb: X = {X:.2f}, Y = {Y:.2f}, Z = {Z:.2f}")

M_A: 3 x 81
jacket under the bulb: X = 35.16, Y = 32.83, Z = 13.01

The jacket's Y = 32.83: it reflects about a third of the light the eye weighs. Its 81 numbers have become three, and that loss is the whole story.

Two spectra, one color

Under the bulb the swatches are one color, but the curves under them have different shapes.

M_A R = XYZ is three equations for 81 unknowns. Any vector b with M_A b = 0 can be added to a reflectance without changing its X, Y, Z: the lamp and the eye see it as no light at all, so colorimetrists call it a metameric black. These vectors form the null space of M_A, with 81 − 3 = 78 dimensions.

The singular value decomposition, np.linalg.svd, finds it. It writes M_A = U Σ Vᵀ, where the 81 rows of Vᵀ are orthonormal directions in reflectance space. M_A stretches the first three by the singular values and sends the other 78 to zero. Those 78 rows span the null space, and as columns they form the matrix N.

The trousers should look like a real cloth, so I start from a second measured textile, a mid green (sample 47), and make the smallest change to it that hides its difference from the lamp. The difference green − jacket splits into a part the lamp sees and a part in the null space. Nᵀ takes the 78 coordinates of the difference along the null-space directions, N puts them back together, and the seen part drops out:

\[R_\text{trousers} = R_\text{jacket} + N N^\mathsf{T}\,(R_\text{green} - R_\text{jacket}).\]

Read the other way, the trousers are the green cloth minus the part of green − jacket that the lamp sees.

U, s, Vt = np.linalg.svd(M_A)
N = Vt[3:].T                                         # 81 x 78, orthonormal columns: what the bulb cannot see
print(f"singular values {s[0]:.2f}, {s[1]:.2f}, {s[2]:.2f}; null space: {N.shape[1]} dimensions")

black = N @ (N.T @ (green - jacket))                 # the metameric black
trousers = jacket + black
print(f"largest |M_A · black|: {np.abs(M_A @ black).max():.0e}")
print(f"trousers' reflectance from {trousers.min():.3f} to {trousers.max():.3f}, "
      f"{np.count_nonzero(np.diff(np.sign(trousers - jacket)))} crossings with the jacket")
for name, R in [("jacket", jacket), ("trousers", trousers), ("green", green)]:
    print(f"{name:9s} X, Y, Z under the bulb: " + ", ".join(f"{v:6.2f}" for v in M_A @ R))

singular values 27.35, 8.51, 8.30; null space: 78 dimensions
largest |M_A · black|: 3e-16
trousers' reflectance from 0.145 to 0.811, 4 crossings with the jacket
jacket    X, Y, Z under the bulb:  35.16,  32.83,  13.01
trousers  X, Y, Z under the bulb:  35.16,  32.83,  13.01
green     X, Y, Z under the bulb:  38.59,  37.40,  10.28

The bulb sees the black as 3e-16, rounding error against values near 30, and jacket and trousers share X, Y, Z to every printed digit while the green cloth is off by 2.7 to 4.6 in each. The trousers' reflectance stays between 0.145 and 0.811, so a fabric can have it. A darker jacket can push it below zero, which a Variation fixes. Two metamers always cross at least three times, and these cross four, circled in the still.

Show code

cross = np.flatnonzero(np.diff(np.sign(trousers - jacket)))
fig, (top, bottom) = plt.subplots(2, 1, figsize=(7, 4.2), sharex=True, gridspec_kw=dict(height_ratios=[1.3, 1]))
top.plot(lam, jacket, color=INK)
top.plot(lam, trousers, color=ACCENT)
for i in cross:                                      # where the difference changes sign, by linear interpolation
    d0, d1 = (trousers - jacket)[i], (trousers - jacket)[i + 1]
    x = lam[i] + 5 * d0 / (d0 - d1)
    top.plot(x, np.interp(x, lam, jacket), "o", ms=6, mfc="none", color=INK)
top.text(450, 0.13, "jacket", color=INK)
top.text(640, 0.72, "trousers", color=ACCENT)
top.set(ylabel="reflectance", ylim=(0, 1))
bottom.axhline(0, color=MUTED, lw=1)
bottom.plot(lam, black, color=ACCENT)
bottom.text(395, 0.14, "metameric black: trousers − jacket", color=ACCENT)
bottom.set(xlabel="λ / nm", ylabel="difference", xlim=(380, 780), ylim=(-0.25, 0.25))
plt.show()

Top: reflectance of the jacket and the trousers from 380 to 780 nm; the curves differ and cross four times, marked with circles. Bottom: their difference, the metameric black, swinging between about −0.17 and +0.21 around zero.

Change the lamp and they part

Switch to daylight and the trousers turn greenish, ΔE00 = 4.98. Under the white LED the difference is 3.49.

The null space belongs to M_A. Daylight builds a different matrix with its own null space, and the black that the bulb could not see is visible to it.

Measuring how different the swatches look needs a distance, and X, Y, Z do not give one: equal steps in them do not look equal. The CIELAB space comes closer. L* is lightness, a* runs from green to red, and b* from blue to yellow, all against the white of the lamp, because the eye adapts to its light: under the warm bulb a gray cloth still looks gray, not orange. ΔE00 (CIEDE2000, the CIE's formula of 2001) is the distance in Lab with corrections for lightness, chroma, and hue, so that equal values look about equally different. About 1 is just noticeable side by side, and above about 2 most people see it.

The swatches in the game and in the still below show each color as the adapted eye sees it, moved from the lamp's white to the screen's daylight white (colour.XYZ_to_sRGB). How the game computes it explains the steps.

def lab(XYZ, lamp):
    """CIELAB against the lamp's own white, the eye adapted to the lamp."""
    white = tristimulus_matrix(lamp).sum(axis=1)     # X, Y, Z of a perfect white reflector
    return colour.XYZ_to_Lab(XYZ / 100, colour.XYZ_to_xy(white / 100))

dE_eye = {}
print(f"{'lamp':7s} {'jacket L*, a*, b*':>20s}   {'trousers L*, a*, b*':>20s}   ΔE00")
for name, lamp in lamps.items():
    M = tristimulus_matrix(lamp)
    lab_j, lab_t = lab(M @ jacket, lamp), lab(M @ trousers, lamp)
    dE_eye[name] = colour.delta_E(lab_j, lab_t, method="CIE 2000")
    row = lambda v: " ".join(f"{x:6.1f}" for x in v)
    print(f"{name:7s} {row(lab_j)}   {row(lab_t)}   {dE_eye[name]:4.2f}")

lamp       jacket L*, a*, b*    trousers L*, a*, b*   ΔE00
A         64.0   -2.9   -5.0     64.0   -2.9   -5.0   0.00
LED-B3    64.0   -3.6   -3.6     63.8   -6.6   -2.7   3.49
D65       64.7   -5.0   -2.5     65.0   -9.5   -0.5   4.98

Under the bulb both have a* = −2.9. In daylight the jacket's a* is −5.0 and the trousers' −9.5, further toward green, while their lightness differs by only 0.4. That green is the black the bulb could not see.

Show code

LAMP_NAMES = {"A": "warm bulb", "LED-B3": "white LED", "D65": "daylight"}

def screen(XYZ, lamp):
    """sRGB as the adapted eye sees it: Bradford from the lamp's white to the sRGB white."""
    white = colour.XYZ_to_xy(tristimulus_matrix(lamp).sum(axis=1) / 100)
    return np.clip(colour.XYZ_to_sRGB(XYZ / 100, illuminant=white, chromatic_adaptation_transform="Bradford"), 0, 1)

fig, (left, right) = plt.subplots(1, 2, figsize=(7, 3.6), sharey=True, gridspec_kw=dict(width_ratios=[1, 1.2]))
for row, (name, lamp) in enumerate(lamps.items()):  # one row per lamp, in both panels
    M = tristimulus_matrix(lamp)
    for col, R in enumerate([jacket, trousers]):
        left.add_patch(plt.Rectangle((col * 1.05, row - 0.35), 1, 0.7, color=screen(M @ R, lamp)))
    left.text(-0.1, row, LAMP_NAMES[name], ha="right", va="center", color=INK)
for col, text in enumerate(["jacket", "trousers"]):
    left.text(col * 1.05 + 0.5, -0.45, text, ha="center", va="bottom", color=INK)
left.set(xlim=(-1.3, 2.1), ylim=(2.6, -0.8))
left.set_axis_off()
right.barh(range(len(lamps)), list(dE_eye.values()), color=ACCENT, height=0.55)
for y, v in enumerate(dE_eye.values()):
    right.text(v + 0.1, y, f"{v:.2f}", va="center", color=ACCENT)
right.axvline(2, color=MUTED, ls="--", lw=1)
right.text(2.08, -0.55, "noticeable", color=MUTED, va="center")
right.set(xlabel="ΔE00, jacket against trousers", xlim=(0, 6))
right.tick_params(axis="y", left=False, labelleft=False)
right.spines["left"].set_visible(False)
right.grid(axis="y", visible=False)
plt.show()

Left: jacket and trousers swatches in screen colors under the warm bulb, the white LED, and daylight, identical under the bulb, the trousers greener in daylight. Right: ΔE00 between them, 0.00 under the bulb, 3.49 under the LED, 4.98 in daylight, against a dashed line at 2, noticeable.

The camera disagrees with the eye

Switch the viewer to camera. In the shop the photo shows a difference at about the threshold, ΔE00 = 2.04, where the eye sees none, and under the LED the camera sees less than the eye, 1.84 against 3.49.

The camera is a Nikon D5100 whose three filters were measured at the UK's National Physical Laboratory (Darrodi et al. 2015). They make a different 3 × 81 matrix with a different null space. To report a color, the camera maps its R, G, B to X, Y, Z with a 3 × 3 matrix, here fitted by least squares on the 99 CIE samples (Least squares). The fit is redone for each lamp with the samples under that lamp, which is the camera's white balance, its counterpart of the eye's adaptation. So the 2.04 in the shop is no mismatch of lamps: the camera's null space is not the eye's.

The 3 × 3 matrix cannot repair that: it only mixes R, G, B after the filters, so a black the filters miss stays invisible, and one they catch stays visible. The camera sees what the eye sees only if its filters are linear combinations of the matching functions, which is called the Luther condition.

def camera_fit(lamp):
    """The 3 x 3 matrix F with X, Y, Z ≈ F · (R, G, B), fitted on the 99 samples under this lamp."""
    k = 100 / (lamp * cmf[:, 1]).sum()
    C = k * (camera * lamp[:, None]).T               # 3 x 81, like M, with the camera's curves
    rgb, xyz = (C @ samples).T, (tristimulus_matrix(lamp) @ samples).T
    F_T, *_ = np.linalg.lstsq(rgb, xyz, rcond=None)  # xyz ≈ rgb @ F_T, one row per sample
    error = colour.delta_E(lab(xyz, lamp), lab(rgb @ F_T, lamp), method="CIE 2000")
    return C, F_T.T, error

F, dE_camera = {}, {}
print(f"{'lamp':7s}  fit error, mean ΔE00   eye ΔE00   camera ΔE00")
for name, lamp in lamps.items():
    C, F[name], error = camera_fit(lamp)
    seen = [lab(F[name] @ C @ R, lamp) for R in (jacket, trousers)]
    dE_camera[name] = colour.delta_E(*seen, method="CIE 2000")
    print(f"{name:7s}  {error.mean():20.2f}   {dE_eye[name]:8.2f}   {dE_camera[name]:11.2f}")

lamp     fit error, mean ΔE00   eye ΔE00   camera ΔE00
A                        1.17       0.00          2.04
LED-B3                   0.88       3.49          1.84
D65                      1.09       4.98          2.94

Even the best matrix misses an average sample by ΔE00 0.88 to 1.17, depending on the lamp, so the D5100's filters fail the Luther condition.

Show code

fig, (left, right) = plt.subplots(1, 2, figsize=(7, 3.6), gridspec_kw=dict(width_ratios=[1, 1.1]))
for j, (band, ls) in enumerate(zip("RGB", ["-", "--", ":"])):
    left.plot(lam, camera[:, j], color=SECOND, ls=ls)
    i = np.argmax(camera[:, j])
    left.text(lam[i], camera[i, j] + 0.04, band, color=SECOND, ha="center")
left.plot(lam, cmf[:, 1], color=INK, ls="--", lw=1.2)
left.text(560, 1.12, "ȳ", color=INK, ha="center")
left.set(xlabel="λ / nm", ylabel="relative sensitivity", xlim=(380, 780), ylim=(0, 1.25))
x = np.arange(len(lamps))
for dx, values, c, who in [(-0.19, dE_eye, ACCENT, "eye"), (0.19, dE_camera, SECOND, "camera")]:
    right.bar(x + dx, list(values.values()), width=0.36, color=c)
    right.text(x[-1] + dx, 0.15, who, color="white", rotation=90, ha="center", va="bottom")  # names inside the daylight bars
    for xi, v in zip(x + dx, values.values()):
        right.text(xi, v + 0.08, f"{v:.2f}", color=c, ha="center", va="bottom",   # white behind: the line at 2
                   bbox=dict(boxstyle="square,pad=0.05", facecolor="white", edgecolor="none"))
right.axhline(2, color=MUTED, ls="--", lw=1, zorder=0)
right.set(xticks=x, xticklabels=["bulb", "LED", "daylight"], ylabel="ΔE00", ylim=(0, 5.8))
right.grid(axis="x", visible=False)
plt.show()

Left: the Nikon D5100's red, green, and blue sensitivities against wavelength, with the matching function ȳ for comparison. Right: ΔE00 between the swatches per lamp for the eye and the camera: warm bulb 0.00 and 2.04, white LED 3.49 and 1.84, daylight 4.98 and 2.94, with a dashed line at 2.

How the game computes it

The state. Two strings, which lamp is on and who looks. {...DEFAULTS} copies the defaults, as dict(DEFAULTS) does in Python:

const DEFAULTS = {lamp: "A", viewer: "eye"};
const state = {...DEFAULTS};

One change. Nothing in the scene moves on its own, so there is no animation loop. A button or a key calls set, which updates the state and calls render, and render computes everything it draws from the state: for each swatch the light S·R it reflects, X, Y, Z, Lab against the lamp's white, and the screen color, then ΔE00 between the two. X, Y, Z are the sums of the spectrum section, one loop over the 81 wavelengths:

function tristimulus(R, S, sens) {
  let norm = 0;
  for (let i = 0; i < S.length; i++) norm += S[i] * DATA.cmf[i][1];
  const k = 100 / norm;                             // a perfect white reflector gets Y = 100
  const out = [0, 0, 0];
  for (let i = 0; i < S.length; i++) {
    for (let j = 0; j < 3; j++) out[j] += k * S[i] * R[i] * sens[i][j];
  }
  return out;
}

For the eye, sens is x̄, ȳ, z̄. For the camera it is the camera's three curves, and the matrix fitted under the current lamp turns its R, G, B into X, Y, Z:

  if (s.viewer === "eye") return tristimulus(R, S, DATA.cmf);
  return mul(DATA.camera_to_xyz[s.lamp], tristimulus(R, S, DATA.camera));

ΔE00. CIEDE2000 takes about 20 lines of JavaScript, most of them about hue angles, which must be subtracted and averaged the short way around the circle. Its last line weighs the three differences, in lightness, chroma, and hue, and a term RT that couples chroma and hue in the blues:

  return Math.sqrt((dL / SL) ** 2 + (dC / SC) ** 2 + (dH / SH) ** 2 + RT * (dC / SC) * (dH / SH));

A line-by-line Python copy checks it against colour.delta_E on every ordered pair of the 99 samples under each lamp and on the game's six settings:

def delta_e00(lab1, lab2):
    """The game's deltaE00, line by line."""
    (L1, a1, b1), (L2, a2, b2) = lab1, lab2
    cos, sin = (lambda d: np.cos(np.radians(d))), (lambda d: np.sin(np.radians(d)))
    hue = lambda b, a: 0.0 if a == 0 and b == 0 else (np.degrees(np.arctan2(b, a)) + 360) % 360
    C7 = ((np.hypot(a1, b1) + np.hypot(a2, b2)) / 2) ** 7
    G = 0.5 * (1 - np.sqrt(C7 / (C7 + 25**7)))
    C1, C2 = np.hypot((1 + G) * a1, b1), np.hypot((1 + G) * a2, b2)
    h1, h2 = hue(b1, (1 + G) * a1), hue(b2, (1 + G) * a2)
    dh = h2 - h1
    if C1 * C2 == 0:
        dh = 0
    elif dh > 180:
        dh -= 360
    elif dh < -180:
        dh += 360
    dL, dC, dH = L2 - L1, C2 - C1, 2 * np.sqrt(C1 * C2) * sin(dh / 2)
    L, C, h = (L1 + L2) / 2, (C1 + C2) / 2, h1 + h2
    if C1 * C2 != 0:
        h = (h + (360 if h < 360 else -360) if abs(h1 - h2) > 180 else h) / 2
    T = 1 - 0.17 * cos(h - 30) + 0.24 * cos(2 * h) + 0.32 * cos(3 * h + 6) - 0.2 * cos(4 * h - 63)
    SL = 1 + 0.015 * (L - 50) ** 2 / np.sqrt(20 + (L - 50) ** 2)
    SC, SH = 1 + 0.045 * C, 1 + 0.015 * C * T
    RT = -2 * np.sqrt(C**7 / (C**7 + 25**7)) * sin(60 * np.exp(-(((h - 275) / 25) ** 2)))
    return np.sqrt((dL / SL) ** 2 + (dC / SC) ** 2 + (dH / SH) ** 2 + RT * (dC / SC) * (dH / SH))

gap, pairs = 0.0, 0
for name, lamp in lamps.items():
    L = lab((tristimulus_matrix(lamp) @ samples).T, lamp)
    for i in range(99):
        for j in range(99):
            if i != j:
                gap = max(gap, abs(delta_e00(L[i], L[j]) - colour.delta_E(L[i], L[j], method="CIE 2000")))
                pairs += 1
    for viewer, M in [("eye", tristimulus_matrix(lamp)), ("camera", F[name] @ camera_fit(lamp)[0])]:
        pair = [lab(M @ R, lamp) for R in (jacket, trousers)]
        gap = max(gap, abs(delta_e00(*pair) - colour.delta_E(*pair, method="CIE 2000")))
print(f"{pairs:,} sample pairs and 6 settings: largest gap to colour.delta_E {gap:.0e}")

29,106 sample pairs and 6 settings: largest gap to colour.delta_E 6e-14

The screen. A screen mixes every color from red, green, and blue pixels, the primaries, whose colors and daylight white the sRGB standard fixes. The Bradford transform comes first: it turns X, Y, Z into three cone-like responses, scales each by the ratio of daylight's white to the lamp's, and turns them back, so a gray under the bulb stays gray on the screen. A second 3 × 3 matrix then says how much of each primary makes those X, Y, Z. The notebook ships their product, one matrix per lamp.

The game applies it, clips to [0, 1], and encodes each channel with the sRGB curve, close to a power of 1/2.4, which spends more of the 256 levels on the dark tones:

    return Math.round(255 * (c <= 0.0031308 ? 12.92 * c : 1.055 * c ** (1 / 2.4) - 0.055));

A color that needs a primary below zero or above full lies outside the screen's gamut, and the game writes "out of screen gamut" under the number. None of the six colors here needs it. A faint tint of the lamp's own color on the wall tells you which lamp is on.

From the notebook into the game. The data travel inside the page, as JSON in a <script type="application/json" id="data"> element: the wavelengths, the two reflectances, the three lamps, x̄, ȳ, z̄, the camera's curves, and the two sets of matrices. Nothing is fetched, so the file runs from disk, and Canvas 2D draws the scene with no library. The cell below fixes the digits of every array and the order of the keys, so the file is byte for byte the same on every run. It also escapes </ as <\/, since a </script> in the data would end the element. The HTML around the game sits in the hidden cell above.

Show code

# The page around the game: the site's light and dark colors, the canvas, the readout, and the buttons.
TEMPLATE = """<!doctype html>
<html lang="en">
<head>
<meta charset="utf-8">
<meta name="viewport" content="width=device-width, initial-scale=1">
<title>Metamers: two gray swatches under three lamps</title>
<style>
:root { --bg: #ffffff; --fg: #1f2a44; --ink: #1f2a44; --accent: #c8553d; --second: #2a7f9e; --muted: #8a8f98; }
@media (prefers-color-scheme: dark) {
  :root { --bg: #141a2a; --fg: #e9ecf3; --ink: #c9d0e3; --accent: #e3826b; --second: #6ec1de; --muted: #8e97ab; }
}
html, body { margin: 0; height: 100%; overflow: hidden; background: var(--bg); color: var(--fg);
  font: clamp(14px, 1.9vw, 22px)/1.25 system-ui, sans-serif; -webkit-text-size-adjust: 100%; }
body { display: flex; flex-direction: column; }
canvas { display: block; width: 100%; flex: none; }
#panel { flex: none; display: flex; flex-wrap: wrap; align-items: center; justify-content: center; gap: 0.4em 1.4em; padding: 0.3em 0.6em 0.6em; }
#readout { color: var(--accent); font-variant-numeric: tabular-nums; white-space: nowrap; }
#de { font-size: 1.8em; font-weight: 600; margin: 0 0.15em; }
#gamut { display: block; color: var(--muted); font-size: max(14px, 0.8em); }
.group { display: flex; align-items: center; gap: 0.35em; }
.group .name { color: var(--muted); margin-right: 0.2em; }
button { font: inherit; line-height: 1.15; color: var(--fg); background: transparent; border: 1.5px solid var(--second);
  border-radius: 4px; padding: 0.25em 0.6em; cursor: pointer; touch-action: manipulation; }
button small { display: block; font-size: max(14px, 0.8em); color: var(--muted); }
button[aria-pressed="true"] { background: var(--second); color: var(--bg); }
button[aria-pressed="true"] small { color: var(--bg); }
@media (pointer: coarse) { button { min-height: 44px; min-width: 44px; } }
@media (max-width: 599px) {
  #panel { flex-direction: column; align-items: stretch; gap: 0.35em; }
  #readout { text-align: center; } .group button { flex: 1; } .group .name { min-width: 3.6em; }
}
</style>
</head>
<body>
<canvas id="scene"></canvas>
<div id="panel">
  <div id="readout">ΔE00<b id="de"></b><span id="verdict"></span><span id="gamut"></span></div>
  <div class="group" role="group" aria-label="Lamp"><span class="name">Lamp</span>
    <button data-lamp="A" title="Warm bulb, key 1">warm bulb<small>2,856 K</small></button>
    <button data-lamp="LED-B3" title="White LED, key 2">white LED<small>4,100 K</small></button>
    <button data-lamp="D65" title="Daylight, key 3">daylight<small>6,500 K</small></button></div>
  <div class="group" role="group" aria-label="Viewer"><span class="name">Viewer</span>
    <button data-viewer="eye" title="Eye, key V">eye</button>
    <button data-viewer="camera" title="Camera, key V">camera</button></div>
</div>
<script type="application/json" id="data">{{DATA}}</script>
<script>{{SCRIPT}}</script>
</body>
</html>
"""

srgb = colour.RGB_COLOURSPACES["sRGB"]

def to_srgb(lamp):
    """Linear sRGB from X, Y, Z / 100, after the Bradford adaptation from the lamp's white to sRGB's."""
    white = colour.xy_to_XYZ(colour.XYZ_to_xy(tristimulus_matrix(lamp).sum(axis=1) / 100))
    cat = colour.adaptation.matrix_chromatic_adaptation_VonKries(
        white, colour.xy_to_XYZ(srgb.whitepoint), transform="Bradford")
    return srgb.matrix_XYZ_to_RGB @ cat

def tint(lamp):
    """The lamp's own color, unadapted, at full brightness, as a hex string."""
    rgb = srgb.matrix_XYZ_to_RGB @ tristimulus_matrix(lamp).sum(axis=1)
    rgb = colour.cctf_encoding(np.clip(rgb / rgb.max(), 0, 1))
    return "#" + "".join(f"{round(255 * c):02x}" for c in rgb)

def rounded(values, digits):
    return np.round(np.asarray(values, dtype=float), digits).tolist()

game_data = {
    "lambda": lam.tolist(),
    "cmf": rounded(cmf, 6),
    "camera": rounded(camera, 5),
    "lamps": {name: rounded(lamp, 4) for name, lamp in lamps.items()},
    "reflectance": {"jacket": rounded(jacket, 5), "trousers": rounded(trousers, 5)},
    "camera_to_xyz": {name: rounded(F[name], 6) for name in lamps},
    "to_srgb": {name: rounded(to_srgb(lamp), 6) for name, lamp in lamps.items()},
    "tint": {name: tint(lamp) for name, lamp in lamps.items()},
}
# the shipped matrices give what colour.XYZ_to_sRGB gives
assert all(np.allclose(colour.cctf_encoding(to_srgb(lamp) @ (tristimulus_matrix(lamp) @ R / 100)),
                       screen(tristimulus_matrix(lamp) @ R, lamp)) for lamp in lamps.values() for R in (jacket, trousers))
data = json.dumps(game_data, sort_keys=True, separators=(",", ":")).replace("</", "<\\/")
script = Path("minigame/metamers/game.js").read_text(encoding="utf-8")
Path("assets").mkdir(exist_ok=True)
Path("assets/metamers.html").write_text(TEMPLATE.replace("{{DATA}}", data).replace("{{SCRIPT}}", script),
                                        encoding="utf-8", newline="\n")
print(f"assets/metamers.html: {Path('assets/metamers.html').stat().st_size / 1000:.0f} kB")

assets/metamers.html: 20 kB

Variations

Match in daylight instead. Take the null space of the daylight matrix in place of M_A. The pair then parts under the bulb, ΔE00 = 3.86, and only barely under the LED, 1.65.
A darker jacket. Make the charcoal sample 74 the jacket and project it toward the gray sample 62: the result dips to −0.046 at 455 nm, a reflectance no cloth can have. minimize(lambda R: np.sum((R - jacket) ** 2), np.clip(projected, 0, 1), method="SLSQP", bounds=[(0, 1)] * 81, constraints={"type": "eq", "fun": lambda R: M_A @ (R - dark)}) finds the closest metamer that stays inside [0, 1] (Minimization with scipy.optimize.minimize). Its minimum is 0.0, and in daylight the pair differs by ΔE00 = 21.7.
Color blindness as two rows. Replace x̄, ȳ, z̄ by the CIE 2006 cone fundamentals and drop the row of the long-wavelength cones, as for a protanope. The null space grows to 79 dimensions, and pairs you can tell apart match for them.
Why a white LED renders red poorly. The color rendering index scores how far test colors shift between a lamp and a glowing black body of the same color temperature, 100 for no shift. Ra averages eight pastel samples, R9 is a saturated red. The white LED scores Ra 84.8 but R9 23.8 (colour.colour_rendering_index), and the bulb, a black body itself, 100. Add that red, SDS_TCS["CIE 1995"]["TCS09"] from colour.quality.cri, as a third swatch and a black body at the LED's 4,102 K as a fourth lamp: the red's chroma C* drops from 70.9 to 62.7 under the LED, ΔE00 = 2.67.

Play: switch the lamp over two gray swatches

A spectrum becomes three numbers

Two spectra, one color

Change the lamp and they part

The camera disagrees with the eye

How the game computes it

Variations

Further reading