Colorimetry

Colorimetry is the measurement of colour as a standard human observer sees it. It reduces a light’s spectrum, a function of wavelength with hundreds of degrees of freedom, to three numbers, the tristimulus values $X$, $Y$ and $Z$, chosen so that two lights with equal tristimulus values look the same to that observer. Everything else in colour measurement is built on these three numbers: chromaticity diagrams, white points, the gamut of a display, the matrices that define RGB colour spaces, and the colour differences used to judge reproduction.

For computer vision the subject matters in two ways. Cameras are not standard observers: their filters see the spectrum differently, so a camera’s raw RGB has to be mapped into a colorimetric space, and two surfaces that a camera separates may look identical to a person, or the reverse. And the standard colour spaces that images are stored in, such as sRGB, are defined colorimetrically, by primaries and a white point in CIE coordinates. Color spaces introduces CIE XYZ in a paragraph and owns sRGB, CIELAB, YCbCr and HSV; radiometry and photometry defines the physical quantities and the luminous weighting that $Y$ shares. This article builds the CIE system itself: colour matching as linear algebra, the standard observers, the tristimulus integral, chromaticity and gamut, illuminants and the Planckian locus, the derivation of an RGB-to-XYZ matrix, and metamerism. Its example computes all of these from an analytic fit to the CIE curves and shows, in a synthetic spectral image, pairs of surfaces that match under one light and break under others.

Definition

Colorimetry assigns to every light, or to every surface under a stated light, three numbers that predict whether two stimuli match in colour for a standard observer. The numbers are linear projections of the light’s spectrum onto three colour-matching functions (CMFs). The CIE’s standard colorimetric observers fix those functions; the CIE’s own definitions are in CIE 015, which was not read for this article. Colorimetry predicts matches, not appearance: it says that two spectra look alike under the same viewing conditions, not what colour either looks like, which also depends on context and adaptation (visionbook, §8.2.1 and §8.3.4).

Intuition

The eye has three kinds of cone, sensitive to long, medium and short wavelengths, and each responds to a weighted sum of the incoming spectrum (visionbook, §8.3.1). Any light, however complicated its spectrum, is therefore reduced to three numbers before vision begins. Two different spectra that produce the same three cone responses are indistinguishable; visionbook (§8.3.3) calls them metameric.

This has a practical consequence that makes colour reproduction possible: a mixture of only three suitable lights can be adjusted to match any colour, and the amounts needed serve as the colour’s coordinates (visionbook, §8.3.2). A display does exactly this with its red, green and blue subpixels. The limitation is that the amounts must be positive: a real light cannot be added in negative quantity, so some colours cannot be matched by any given three primaries. They lie outside those primaries’ gamut (visionbook, §8.3.5).

The CIE system sidesteps the limitation by measuring rather than reproducing. Its three functions are chosen to be non-negative at every wavelength, so every real light gets non-negative coordinates and every device’s primaries can be described in them; the price is that no set of real, non-negative primaries matches these functions (visionbook, §8.3.6).

Formal Definition

Colour matching as linear algebra

Sample the visible range at $N$ wavelengths, so that a light is a vector $\mathbf{t} \in \mathbb{R}^N$. Visionbook (§8.3.1, Eq. 8.1) writes the cone responses as $\mathbf{C}_{\text{eye}}\mathbf{t}$, with the three cone sensitivities in the rows of the $3 \times N$ matrix $\mathbf{C}_{\text{eye}}$. Any three sensors with sensitivities $\mathbf{C} = \mathbf{R}\,\mathbf{C}_{\text{eye}}$, for a full-rank $3 \times 3$ matrix $\mathbf{R}$, measure the same information (Eq. 8.2): their readings determine the cone responses and the reverse. A colour reproduction system reads $\mathbf{C}\mathbf{t}$, mixes its three primaries (the columns of $\mathbf{P}$) with weights $\mathbf{M}\mathbf{C}\mathbf{t}$, and emits $\mathbf{P}\mathbf{M}\mathbf{C}\mathbf{t}$. For the output to match every input in the eye, visionbook derives (Eqs. 8.3–8.6)

$$\mathbf{C} = \mathbf{R}\,\mathbf{C}_{\text{eye}}, \qquad \mathbf{M} = (\mathbf{C}\mathbf{P})^{-1}.$$

Two consequences follow, both in our derivation from these equations. First, the sensitivities that match a given set of primaries, with $\mathbf{M} = \mathbf{I}$, are $\mathbf{C} = (\mathbf{C}_{\text{eye}}\mathbf{P})^{-1}\mathbf{C}_{\text{eye}}$: they are the colour-matching functions of those primaries, the amounts of each primary needed to match a unit of each wavelength. Second, a monochromatic light at a wavelength where any of these functions is negative needs a negative amount of some primary, which is the gamut limit of the intuition above in algebraic form. The example shows that for the sRGB primaries this happens at every sampled wavelength from 380 to 780 nm, computed with the fit described below.

The CIE standard observers

The CIE standardized colour-matching functions derived from human matching experiments of the 1920s (Wyman, Sloan & Shirley, 2013, §1). There are two standard observers, the 1931 observer for a 2° field of view and the 1964 observer for a 10° field (Wyman et al., §1.1); Wyman et al. call the 1931 curves the standard most commonly used in graphics, and sRGB’s white point is given for the 1931 observer (colour-science, srgb.py, which takes D65 from its 2° table). The 1931 functions are written $\bar{x}(\lambda)$, $\bar{y}(\lambda)$ and $\bar{z}(\lambda)$.

Visionbook (§8.3.6) states the design choice that defines XYZ: the CIE functions were made non-negative at every wavelength, and as a result no set of non-negative primaries forms a colour-matching system with them. XYZ is a valid space in which to measure colours, but its “primaries” are not physical lights. Color spaces adds that $Y$ is defined to be proportional to luminance, which is why $Y$ carries the same luminous weighting that radiometry and photometry uses for lumens and candelas.

The 1931 functions are published as tables (CIE 015, not read here). Wyman, Sloan and Shirley (2013) give analytic fits that avoid the tables. Their most accurate one, a sum of piecewise Gaussians (§2.2, Eq. 4),

$$\bar{x}(\lambda) \approx \sum_i \alpha_i \exp\!\Big(-\tfrac{1}{2}\big[(\lambda - \beta_i)\, S(\lambda - \beta_i, \gamma_i, \delta_i)\big]^2\Big),$$

uses three lobes for $\bar{x}$ and two each for $\bar{y}$ and $\bar{z}$, where the selector $S(x, \gamma, \delta)$ is $\gamma$ for $x < 0$ and $\delta$ otherwise, so $\gamma_i$ and $\delta_i$ are inverse widths on the two sides of the centre $\beta_i$ (Table 1 gives the 28 coefficients, Listing 1 the same fit as C code). Measured against the CIE curves sampled at 1 nm, the fit’s largest squared errors are $2.0 \times 10^{-4}$, $6.4 \times 10^{-5}$ and $4.9 \times 10^{-4}$ for $\bar{x}$, $\bar{y}$ and $\bar{z}$, and its mean squared errors $3.1 \times 10^{-5}$, $7.1 \times 10^{-6}$ and $1.6 \times 10^{-5}$ (Table 2, read from the page image). The authors report that this is less than the variation between repeated measurements on a single subject in the experiments behind the 1964 observer, and similar to linearly interpolating 10 nm samples of the curves; they also note that the CIE curves themselves were interpolated to 1 nm (§3). The example uses this fit throughout, so its numbers carry the fit’s error, which it measures where an exact value is known.

Tristimulus values

For a surface of spectral reflectance $R(\lambda)$ under a light of spectral power $S(\lambda)$, the tristimulus values are

$$X = k\!\int S(\lambda) R(\lambda)\, \bar{x}(\lambda)\, d\lambda, \quad Y = k\!\int S(\lambda) R(\lambda)\, \bar{y}(\lambda)\, d\lambda, \quad Z = k\!\int S(\lambda) R(\lambda)\, \bar{z}(\lambda)\, d\lambda,$$

with $R = 1$ for a light viewed directly. Visionbook (§8.2.3) models diffuse reflection as the wavelength-by-wavelength product of the incident spectrum and the surface reflectance, which is the $S R$ inside the integrals. The constant $k$ sets the scale. colour-science’s sd_to_XYZ_integration docstring (v0.4.7) gives the two conventions: for reflecting or transmitting objects, $k$ is chosen so that $Y = 100$ for a perfect reflector ($R = 1$ everywhere), that is $k = 100/\!\int S\bar{y}\,d\lambda$; for self-luminous sources $k$ is chosen for convenience, and when $Y$ is to equal a photometric quantity numerically, $k$ is the maximum luminous efficacy $K_m = 683$ lm/W ($Y$ is then a luminance for spectral radiance), as radiometry and photometry explains. Relative colorimetry, with $Y = 100$ or $Y = 1$ for white, is the convention of colour spaces and of this article.

The integral is computed on sampled data. colour-science offers plain integration and the ASTM E308 method, which uses weighting tables; its sd_to_XYZ defaults to “ASTM E308”, so a cross-check against it must name the method. The example uses the trapezoid rule at 1 nm.

Chromaticity and the chromaticity diagram

The overall level of a colour is carried by the scale of $(X, Y, Z)$; its chromaticity is what remains after normalizing it away (visionbook, §8.3.6; color spaces):

$$x = \frac{X}{X + Y + Z}, \qquad y = \frac{Y}{X + Y + Z}, \qquad z = 1 - x - y .$$

A colour can be given as $(x, y, Y)$, chromaticity plus luminance. Mixtures behave simply in this plane (our derivation): adding two lights adds their tristimulus values, and the chromaticity of the sum is the average of the two chromaticities weighted by $X + Y + Z$ of each light. A mixture of two lights therefore lies on the straight segment between their chromaticities, and every mixture of three lights with non-negative weights lies in the triangle they span.

Monochromatic lights trace the spectral locus, the horseshoe-shaped curve of the chromaticity diagram. Every physical light is a non-negative mixture of monochromatic lights, so all physical chromaticities lie in the convex hull of the locus, the region bounded by the locus and the straight segment joining its two ends, which holds mixtures of the shortest and longest wavelengths and contains no monochromatic light (our derivation from the mixing rule). The chromaticity diagram is a projection, not a perceptual map: equal distances in it are not equal visible differences, which is why CIELAB and its colour differences exist (color spaces).

Primaries, gamut and RGB matrices

A display or a colour space with three primaries of chromaticities $(x_r, y_r)$, $(x_g, y_g)$, $(x_b, y_b)$ can produce exactly the chromaticities inside their triangle; that triangle is its gamut in the chromaticity diagram. The spectral locus is curved outwards, so a triangle with corners inside it leaves parts of the locus outside, and saturated colours near the locus cannot be produced (our derivation). sRGB’s primaries are $(0.640, 0.330)$, $(0.300, 0.600)$ and $(0.150, 0.060)$ with the D65 white point (colour-science srgb.py, citing IEC 61966-2-1, which was not read; color spaces gives the same values).

A linear RGB space is defined colorimetrically by a $3 \times 3$ matrix that maps linear $(R, G, B)$ to $(X, Y, Z)$. The primaries fix its columns up to scale; the white point fixes the scales. colour-science computes this normalised primary matrix (NPM) in normalised_primary_matrix, citing SMPTE RP 177 (not read): put each primary’s $(x, y, 1 - x - y)$ in the columns of $\mathbf{P}$, convert the white point to $\mathbf{W} = (x_w/y_w,\, 1,\, (1 - x_w - y_w)/y_w)$, solve $\mathbf{P}\mathbf{s} = \mathbf{W}$ for the column scales, and form

$$\mathbf{M}_{\text{RGB}\to\text{XYZ}} = \mathbf{P}\,\operatorname{diag}(\mathbf{s}),$$

so that $R = G = B = 1$ gives the white with $Y = 1$. Its middle row gives luminance as a weighted sum of linear $R$, $G$ and $B$; for sRGB these are the weights 0.2126, 0.7152 and 0.0722 that color spaces prints. The inverse matrix takes XYZ to linear RGB, and a colour outside the gamut comes out with a negative component, or above 1 if it is too bright. A camera’s colour correction matrix plays the same role for the camera’s own RGB, in our reading fitted rather than derived: the camera’s sensitivities do not match the CIE functions (image formation), and unless they are exact linear combinations of them, no matrix maps every camera colour exactly.

White points and illuminants

Colorimetry needs a reference white. For surfaces it is the perfect reflector under the light in use, which is why the light must be stated; for colour spaces it is a chromaticity written into the definition. colour-science’s table of illuminant chromaticities (v0.4.7, citing CIE 015) gives, for the 1931 2° observer:

• E, the equal-energy spectrum, constant in wavelength: exactly $(1/3, 1/3)$, which means (our derivation) that the three 1931 functions have equal areas.
• A, at $(0.44758, 0.40745)$, described in colour-science’s sd_CIE_standard_illuminant_A as typical domestic tungsten-filament lighting and as a Planckian radiator at approximately 2856 K.
• D65, at $(0.31270, 0.32900)$, the white of sRGB.

The same table gives different values for the 1964 10° observer, for example $(0.31382, 0.33100)$ for D65 and $(0.45117, 0.40594)$ for A: a white point belongs to an observer as well as to a light.

Planckian light and illuminant A

A blackbody’s spectrum depends only on its temperature (Planck’s law, stated in radiometry and photometry), so its chromaticity traces a curve in the diagram as the temperature varies, the Planckian locus; the example prints seven points on it, from $(0.581, 0.398)$ at 1500 K to $(0.257, 0.258)$ at 20000 K. With the second radiation constant $c_2 = hc/k$, the spectral shape is

$$B_\lambda(T) \propto \frac{\lambda^{-5}}{\exp\!\big(c_2/(\lambda T)\big) - 1},$$

and it depends on $c_2$ and $T$ only through their ratio. That detail explains a standard’s temperature. colour-science’s planck_law docstring says that colorimetry should use the $c_2$ of the current International Temperature Scale, ITS-90, $1.4388 \times 10^{-2}$ m K. Illuminant A, however, is defined by a formula with an older constant: colour-science computes it, citing CIE 015:2004, as

$$S_A(\lambda) = 100 \left(\frac{560}{\lambda}\right)^{5} \frac{\exp\!\big(1.435 \times 10^{7}/(2848 \times 560)\big) - 1}{\exp\!\big(1.435 \times 10^{7}/(2848\, \lambda)\big) - 1}, \quad \lambda \text{ in nm},$$

a Planck spectrum with $c_2 = 1.435 \times 10^{-2}$ m K at 2848 K, normalized to 100 at 560 nm. Since only $c_2/T$ matters, this is exactly the ITS-90 Planck spectrum at $2848 \times 1.4388/1.435 = 2855.54$ K (our derivation, confirmed to $4 \times 10^{-15}$ by check C3 in the example), which is why A is described as approximately 2856 K.

Metamerism

Two spectra are metamers when they give the same tristimulus values (visionbook, §8.3.3): any two spectra with the same projection onto the colour-matching functions look the same, and since the spectra live in a high-dimensional space and vision sees only a 3-D projection, the space of metamers is large. For surfaces the projection includes the light. Writing $\mathbf{A}_S$ for the $3 \times N$ matrix whose rows are $k\,S(\lambda)\bar{x}(\lambda)$, $k\,S(\lambda)\bar{y}(\lambda)$ and $k\,S(\lambda)\bar{z}(\lambda)$ times the integration weights, the tristimulus values of reflectance $\mathbf{r}$ are $\mathbf{A}_S\mathbf{r}$, and (our derivation)

$$\mathbf{A}_S(\mathbf{r} + \mathbf{b}) = \mathbf{A}_S\mathbf{r} \quad \text{for every } \mathbf{b} \in \operatorname{null}(\mathbf{A}_S).$$

Any reflectance plus any component in the null space of $\mathbf{A}_S$, kept within $[0, 1]$, is a metamer of the original under light $S$. The null space depends on $S$: a component invisible under one light is in general visible under another, so the pair matches under one light and not under the other. This is metamerism with respect to the illuminant. In our reading the same argument applies to the observer: replacing the CMFs, by another person’s or by a camera’s sensitivities, changes the matrix and so the null space, which is how a camera can separate two surfaces that look identical to a person, or merge two that do not.

Properties

• Linearity. Tristimulus values are linear in the spectrum, so they add for mixed lights, and any linear change of basis, such as XYZ to linear RGB, is a $3 \times 3$ matrix (visionbook, §8.3.4). Chromaticity is not linear, but mixtures lie on segments and in triangles.
• Positivity of XYZ. The CIE functions are non-negative, so every physical light has non-negative $X$, $Y$, $Z$; the price is that XYZ has no physical primaries (visionbook, §8.3.6).
• No real triangle covers the locus. Three physical primaries lie inside the spectral locus, and so do their mixtures (our derivation); colour-matching functions of real primaries are negative somewhere.
• A white point needs a light and an observer. E, A and D65 have different chromaticities, and each changes between the 2° and 10° observers.
• Metamerism depends on the light and the observer. A match is a statement about one projection; it can break when the light, the observer or the sensor changes.
• Colorimetry is not appearance. Equal tristimulus values predict a match side by side under the same conditions; the colour seen also depends on adaptation and surroundings (visionbook, §8.3.4), and that is the territory of color constancy and chromatic adaptation.

Examples

Metamers you can build, and the lights that break them

The script below computes the CIE system from Wyman et al.’s multi-lobe fit to the 1931 CMFs, on 380–780 nm at 1 nm steps with the trapezoid rule, with nothing downloaded: no CIE tables, illuminant tables or images. It has five parts.

• M1, the fit and the equal-energy white. The fit is coded twice, once from Table 1 (vectorized) and once by transcribing Listing 1 line by line; check C1 compares them. The chromaticity of the equal-energy spectrum E is then compared with its exact value (1/3, 1/3).
• M2, the Planckian locus and illuminant A. Chromaticities of ITS-90 Planck spectra from 1500 to 20000 K; illuminant A from colour-science’s formula, checked against the values printed in its docstring (C2) and against the ITS-90 Planck spectrum at 2855.54 K (C3); A against scikit-image’s A white and colour-science’s tabulated A.
• M3, the sRGB matrix. The NPM derived from the sRGB primaries and D65 (C4: it maps $R = G = B = 1$ to the white; C4b: its luminance row against Color Spaces’ weights), compared with scikit-image’s xyz_from_rgb and with colour-science’s 4-decimal MATRIX_sRGB_TO_XYZ, together with the D65 whites stored in scikit-image.
• M4, the spectral locus and the sRGB triangle. Each monochromatic stimulus is converted to linear sRGB; it lies inside the triangle only if all three components are non-negative.
• M5, a spectral image of metamers. Six base reflectances (five Gaussian bumps of reflectance 0.2–0.8 centred at 430–630 nm, and a flat grey 0.5) are given null-space components for light L1, a 6500 K blackbody (a daylight-like light, not D65, whose spectrum colour-science builds from tabulated basis functions, data that are not downloaded here). The components come from an 8-term cosine basis: the 5-dimensional null space of the $3 \times 8$ matrix $\mathbf{A}_{S_1}\mathbf{B}$ gives smooth spectra that L1 cannot see. Row $i$ adds the $i$-th null-space direction (the first again for the sixth row), scaled to a maximum of 1 and multiplied by $t\,a_{\max}$, $t = 1/3, 2/3, 1$, where $a_{\max}$ is the largest amplitude that keeps the reflectance in $[0.02, 0.98]$. The 24 reflectances fill a $192 \times 128$ spectral image of $32 \times 32$-pixel patches (rows: base reflectances; columns: base, then the three amplitudes). The image is rendered to XYZ under L1, under CIE A (L2), and under L3, a three-band light made of Gaussians at 450, 540 and 610 nm with σ = 10 nm and equal peaks (our construction, loosely like a white LED made of three narrow emitters). Each light is normalized so that its perfect reflector has $Y = 100$. Each metamer is compared with its base by $\Delta E^*_{ab}$, the Euclidean distance in CIELAB (color spaces), with each light’s own perfect white as the reference, and by whether the two patches get identical 8-bit sRGB codes when the image is rendered with the derived matrix, clipped to $[0, 1]$ and encoded with the sRGB transfer function. No chromatic adaptation is applied, so the chart under A is rendered as orange as the light.

The checks (C1–C6) and data checks (D1–D3) verify the code and the data; the eleven hypotheses (H1–H11) were fixed before any part of the example was run:

• Checks: C1, the two codings of the fit agree to $10^{-12}$; C2, the A formula reproduces colour-science’s docstring values to $10^{-6}$; C3, A equals the ITS-90 Planck shape at 2855.54 K to $10^{-12}$; C4, the NPM maps $(1, 1, 1)$ to the D65 white to $10^{-12}$, and C4b, its luminance row rounds to (0.2126, 0.7152, 0.0722); C5, the metamers’ XYZ under L1 differ from their bases’ by less than $10^{-9}$; C6, our CIELAB agrees with scikit-image’s xyz2lab to $10^{-10}$ (scikit-image divides by the white of the illuminant it is given, so input already divided by our white, with illuminant "E", whose white is $(1, 1, 1)$, gives CIELAB relative to our white).
• Data checks: D1, every reflectance in the image lies in $[0.02, 0.98]$; D2, every patch is uniform and the 24 patches are distinct; D3, the null space has dimension 5 and every row’s $a_{\max}$ exceeds 0.05.
• Hypotheses: (H1) the fit puts E within 0.01 of $(1/3, 1/3)$ in $x$ and $y$; (H2) the fit puts A within 0.005 of colour-science’s and scikit-image’s A in each coordinate; (H3) the ITS-90 Planck spectrum at 2856 K and the A formula differ by less than $10^{-4}$ in $x$ and $y$; (H4) $x$ decreases strictly along the locus from 1500 to 20000 K; (H5a) the NPM differs from scikit-image’s matrix by less than $5 \times 10^{-4}$ in every entry, and (H5b) rounded to 4 decimals equals colour-science’s matrix; (H6) the $Z$ row sum of scikit-image’s matrix differs from scikit-image’s own D65 $Z$ (1.08883) by more than $10^{-5}$; (H7) no monochromatic stimulus between 380 and 780 nm lies inside the sRGB triangle; (H8) every pair’s $\Delta E^*_{ab}$ under L1 is below $10^{-6}$; (H9) under L2 and under L3, $\Delta E^*_{ab}$ increases strictly with $t$ in every row; (H10) the median $\Delta E^*_{ab}$ over the 18 pairs is larger under the three-band light L3 than under A; (H11) under L1 all 18 pairs have identical 8-bit codes, and under L2 and L3 none of the six pairs at $t = 1$ does.

The script runs in about a second with NumPy 2.3.5, SciPy 1.16.3 and scikit-image 0.25.2 (Python 3.13).

import os
os.environ["OPENBLAS_NUM_THREADS"] = "1"   # one BLAS thread, before NumPy
os.environ["OMP_NUM_THREADS"] = "1"
import math
import numpy as np
from scipy.linalg import null_space
from skimage.color import xyz2lab, xyz_tristimulus_values
from skimage.color.colorconv import xyz_from_rgb

lam = np.arange(380.0, 781.0)                  # wavelength grid, nm
w = np.ones_like(lam); w[0] = w[-1] = 0.5      # trapezoid weights, 1 nm steps

# Wyman, Sloan & Shirley (2013), Eq. 4 and Table 1: multi-lobe fit to the CIE 1931 CMFs
AL = {"x": [0.362, 1.056, -0.065], "y": [0.821, 0.286], "z": [1.217, 0.681]}
BE = {"x": [442.0, 599.8, 501.1], "y": [568.8, 530.9], "z": [437.0, 459.0]}
GA = {"x": [0.0624, 0.0264, 0.0490], "y": [0.0213, 0.0613], "z": [0.0845, 0.0385]}
DE = {"x": [0.0374, 0.0323, 0.0382], "y": [0.0247, 0.0322], "z": [0.0278, 0.0725]}

def cmf(l):
    out = []
    for c in "xyz":
        s = np.zeros_like(l)
        for a, b, g, d in zip(AL[c], BE[c], GA[c], DE[c]):
            t = (l - b) * np.where(l < b, g, d)    # gamma below the centre, delta above
            s += a * np.exp(-0.5 * t * t)
        out.append(s)
    return np.stack(out)                       # 3 x N: xbar, ybar, zbar

def listing(wave):                             # Listing 1 of the paper, transcribed line by line
    t1 = (wave-442.0)*(0.0624 if wave < 442.0 else 0.0374)
    t2 = (wave-599.8)*(0.0264 if wave < 599.8 else 0.0323)
    t3 = (wave-501.1)*(0.0490 if wave < 501.1 else 0.0382)
    x = 0.362*math.exp(-0.5*t1*t1) + 1.056*math.exp(-0.5*t2*t2) - 0.065*math.exp(-0.5*t3*t3)
    t1 = (wave-568.8)*(0.0213 if wave < 568.8 else 0.0247)
    t2 = (wave-530.9)*(0.0613 if wave < 530.9 else 0.0322)
    y = 0.821*math.exp(-0.5*t1*t1) + 0.286*math.exp(-0.5*t2*t2)
    t1 = (wave-437.0)*(0.0845 if wave < 437.0 else 0.0278)
    t2 = (wave-459.0)*(0.0385 if wave < 459.0 else 0.0725)
    z = 1.217*math.exp(-0.5*t1*t1) + 0.681*math.exp(-0.5*t2*t2)
    return x, y, z

CMF = cmf(lam)

def weights(S):                                # X = k sum(w S R xbar), k: Y = 100 for R = 1
    k = 100 / np.sum(w * S * CMF[1])
    return k * CMF * (w * S)

def xy(XYZ):
    return XYZ[:2] / XYZ.sum()

def planck(l_nm, T, c2=1.4388e-2):             # spectral shape only; c1 cancels in xy
    l = l_nm * 1e-9
    return l ** -5 / np.expm1(c2 / (l * T))

# CIE standard illuminant A as in colour-science v0.4.7, colour/colorimetry/illuminants.py
# (sd_CIE_standard_illuminant_A), ported. Copyright 2013 Colour Developers.
# Redistribution and use in source and binary forms, with or without modification, are
# permitted provided that the following conditions are met:
# 1. Redistributions of source code must retain the above copyright notice, this list of
# conditions and the following disclaimer.
# 2. Redistributions in binary form must reproduce the above copyright notice, this list
# of conditions and the following disclaimer in the documentation and/or other materials
# provided with the distribution.
# 3. Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific prior
# written permission.
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
# EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
# MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
# THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
# OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
# LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
def cie_A(l):
    return 100 * (560 / l) ** 5 * (np.expm1(1.435e7 / (2848 * 560)) / np.expm1(1.435e7 / (2848 * l)))

# Normalised primary matrix (same computation as colour-science's normalised_primary_matrix,
# same licence as above): columns are the primaries' XYZ, scaled so that RGB = 1 gives the white.
def npm(prim, wxy):
    P = np.vstack([prim[:, 0], prim[:, 1], 1 - prim.sum(1)])
    W = np.array([wxy[0] / wxy[1], 1.0, (1 - wxy.sum()) / wxy[1]])
    return P @ np.diag(np.linalg.solve(P, W)), W

# ---- M1: the CMF port and the equal-energy white ----
print("C1 max|Table 1 fit - Listing 1| over 380-780 nm: %.3e (tol 1e-12)"
      % np.abs(CMF - np.array([listing(v) for v in lam]).T).max())
e = xy(weights(np.ones_like(lam)).sum(1))
print("M1 E white: x %.6f y %.6f  |dx| %.3e |dy| %.3e (H1 bar 0.01)" % (e[0], e[1], abs(e[0]-1/3), abs(e[1]-1/3)))

# ---- M2: Planckian locus and illuminant A ----
doc = np.array([14.7080384, 59.8610989, 100.0, 129.0427389, 198.2612232])
print("C2 max|A formula - colour docstring| at 400,500,560,600,700 nm: %.3e (tol 1e-6)"
      % np.abs(cie_A(np.array([400., 500., 560., 600., 700.])) - doc).max())
Teq = 2848 * 1.4388 / 1.435
pa = planck(lam, Teq) / planck(np.array([560.0]), Teq)
print("C3 max|A / Planck(%.2f K, ITS-90 c2) - 1|: %.3e (tol 1e-12)" % (Teq, np.abs(cie_A(lam) / 100 / pa - 1).max()))
temps = [1500, 2000, 2856, 4000, 6500, 10000, 20000]
loc = np.array([xy(weights(planck(lam, T)).sum(1)) for T in temps])
for T, (x_, y_) in zip(temps, loc):
    print("M2 Planck %5d K: x %.6f y %.6f" % (T, x_, y_))
print("H4 x strictly decreasing along the locus:", bool(np.all(np.diff(loc[:, 0]) < 0)))
a_f = xy(weights(cie_A(lam)).sum(1))
p56 = loc[temps.index(2856)]
XA = xyz_tristimulus_values(illuminant="A", observer="2")
sk_a = XA[:2] / XA.sum()
cs_a = np.array([0.44758, 0.40745])
print("M2 A formula      : x %.6f y %.6f" % tuple(a_f))
print("M2 skimage A      : x %.6f y %.6f   colour-science A: x %.6f y %.6f" % (*sk_a, *cs_a))
print("H2 |A formula - colour-science A| %.3e %.3e; - skimage A %.3e %.3e (bar 0.005)"
      % (*np.abs(a_f - cs_a), *np.abs(a_f - sk_a)))
print("H3 |Planck 2856 K - A formula| %.3e %.3e (bar 1e-4)" % tuple(np.abs(p56 - a_f)))

# ---- M3: the sRGB matrix ----
prim = np.array([[0.64, 0.33], [0.30, 0.60], [0.15, 0.06]])     # sRGB primaries (colour-science srgb.py)
M, W = npm(prim, np.array([0.3127, 0.3290]))                     # D65, 2 degree observer
cs_M = np.array([[0.4124, 0.3576, 0.1805], [0.2126, 0.7152, 0.0722], [0.0193, 0.1192, 0.9505]])
print("M3 derived NPM:\n", np.array2string(M, precision=6, suppress_small=True))
print("C4 max|NPM (1,1,1) - D65 XYZ|: %.3e (tol 1e-12); C4b Y row x1e4:" % np.abs(M.sum(1) - W).max(),
      np.rint(M[1] * 1e4).astype(int).tolist(), "(color spaces: [2126, 7152, 722])")
print("H5a max|NPM - skimage xyz_from_rgb|: %.3e (bar 5e-4)" % np.abs(M - xyz_from_rgb).max())
print("H5b NPM rounded to 4 decimals equals colour-science matrix:",
      bool(np.array_equal(np.rint(M * 1e4).astype(int), np.rint(cs_M * 1e4).astype(int))))
rs = xyz_from_rgb.sum(1)
sk_w = xyz_tristimulus_values(illuminant="D65", observer="2")
print("M3 whites (X, Y, Z): from xy %.6f %.6f %.6f | skimage D65 %.6f %.6f %.6f | skimage matrix row sums %.6f %.6f %.6f"
      % (*W, *sk_w, *rs))
print("H6 |skimage row-sum Z - skimage D65 Z|: %.3e (bar > 1e-5)" % abs(rs[2] - sk_w[2]))

# ---- M4: the spectral locus against the sRGB triangle ----
rgb = np.linalg.solve(M, CMF)                  # linear sRGB of each monochromatic stimulus
inside = np.all(rgb >= 0, axis=0)
print("H7 monochromatic stimuli inside the sRGB triangle: %d of %d; largest min(r,g,b)/(r+g+b): %.4e"
      % (inside.sum(), lam.size, (rgb.min(0) / rgb.sum(0)).max()))
for l0 in (450, 520, 580, 650):
    c = CMF[:, int(l0 - 380)]
    print("M4 %d nm: x %.6f y %.6f" % (l0, *(c[:2] / c.sum())))

# ---- M5: a spectral image of metamer pairs under three lights ----
S1 = planck(lam, 6500.0)                                         # match light (not D65)
S2 = cie_A(lam)                                                  # CIE A
S3 = sum(np.exp(-0.5 * ((lam - c) / 10.0) ** 2) for c in (450.0, 540.0, 610.0))   # three bands
lights = {"L1 6500 K": S1, "L2 CIE A": S2, "L3 3-band": S3}
u = (lam - 380) / 400
B = np.stack([np.cos(np.pi * k * u) for k in range(8)], 1)      # smooth basis, 401 x 8
N = null_space(weights(S1) @ B)                                  # coefficients invisible under L1
base = [0.2 + 0.6 * np.exp(-0.5 * ((lam - c) / 50.0) ** 2) for c in (430., 480., 530., 580., 630.)]
base.append(np.full_like(lam, 0.5))
ts = (1/3, 2/3, 1.0)
refl, amax = [], []
for i, R in enumerate(base):
    b = B @ N[:, i % 5]; b = b / np.abs(b).max()                 # a null-space component: invisible under L1
    up = np.where(b > 0, (0.98 - R) / np.where(b > 0, b, 1), np.inf)
    dn = np.where(b < 0, (R - 0.02) / np.where(b < 0, -b, 1), np.inf)
    a = min(up.min(), dn.min()); amax.append(a)
    refl.append([R] + [R + t * a * b for t in ts])
PS = 32
img = np.zeros((6 * PS, 4 * PS, lam.size))                       # rows: bases; columns: base, t = 1/3, 2/3, 1
for i in range(6):
    for j in range(4):
        img[i*PS:(i+1)*PS, j*PS:(j+1)*PS] = refl[i][j]
patches = img.reshape(6, PS, 4, PS, -1).transpose(0, 2, 1, 3, 4).reshape(24, PS * PS, -1)
print("D1 reflectance min %.6f max %.6f within [0.02, 0.98]:" % (img.min(), img.max()),
      bool(img.min() >= 0.02 - 1e-12 and img.max() <= 0.98 + 1e-12))
print("D2 every patch uniform:", bool(all(np.all(p == p[0]) for p in patches)),
      "; distinct patches:", len({tuple(np.round(p[0], 12)) for p in patches}))
print("D3 null-space dimension %d; a_max per row:" % N.shape[1], " ".join("%.4f" % a for a in amax),
      "all > 0.05:", bool(min(amax) > 0.05))

def lab(XYZ, white):                           # CIELAB with scikit-image's constants
    r = XYZ / white
    f = np.where(r > 0.008856, np.cbrt(r), 7.787 * r + 16 / 116)
    return np.stack([116 * f[..., 1] - 16, 500 * (f[..., 0] - f[..., 1]), 200 * (f[..., 1] - f[..., 2])], -1)

def encode(c):                                 # sRGB transfer function (see color spaces)
    return np.where(c <= 0.0031308, 12.92 * c, 1.055 * np.clip(c, 0, None) ** (1 / 2.4) - 0.055)

maxd, med = 0.0, {}
for name, S in lights.items():
    Wt = weights(S)
    xyz_img = np.tensordot(img, Wt, axes=([2], [1]))             # 192 x 128 x 3 image of XYZ
    pv = xyz_img.reshape(6, PS, 4, PS, 3).mean((1, 3))           # patch values, 6 x 4 x 3
    white = Wt.sum(1)
    L = lab(pv, white)
    maxd = max(maxd, np.abs(L - xyz2lab((pv / white).reshape(-1, 1, 3), illuminant="E").reshape(6, 4, 3)).max())
    dE = np.linalg.norm(L[:, 1:] - L[:, :1], axis=-1)            # 6 x 3: each metamer against its base
    lin = np.linalg.solve(M, (xyz_img / 100).reshape(-1, 3).T).T
    clipped = int(np.sum((lin < 0) | (lin > 1)))
    code = np.rint(255 * encode(np.clip(lin, 0, 1))).astype(int).reshape(6, PS, 4, PS, 3)[:, 0, :, 0]
    same = np.all(code[:, 1:] == code[:, :1], axis=-1)           # 6 x 3
    med[name] = np.median(dE)
    print("%s | white x %.4f y %.4f | max|dXYZ| %.3e | max dE %.3e | median dE %.6f | identical 8-bit pairs %d of 18 (t = 1: %d of 6) | clipped values %d"
          % (name, *xy(white), np.abs(pv[:, 1:] - pv[:, :1]).max(), dE.max(), med[name], same.sum(), same[:, 2].sum(), clipped))
    for i in range(6):
        print("   row %d dE at t = 1/3, 2/3, 1: %.6f %.6f %.6f  increasing: %s"
              % (i, *dE[i], bool(dE[i, 0] < dE[i, 1] < dE[i, 2])))
print("C5 is the L1 max|dXYZ| above (tol 1e-9); C6 max|own Lab - skimage xyz2lab|: %.3e (tol 1e-10)" % maxd)
print("H10 median dE: L3 %.6f > L2 %.6f: %s" % (med["L3 3-band"], med["L2 CIE A"], bool(med["L3 3-band"] > med["L2 CIE A"])))

Output:

C1 max|Table 1 fit - Listing 1| over 380-780 nm: 0.000e+00 (tol 1e-12)
M1 E white: x 0.332998 y 0.333716  |dx| 3.356e-04 |dy| 3.824e-04 (H1 bar 0.01)
C2 max|A formula - colour docstring| at 400,500,560,600,700 nm: 5.539e-08 (tol 1e-6)
C3 max|A / Planck(2855.54 K, ITS-90 c2) - 1|: 3.997e-15 (tol 1e-12)
M2 Planck  1500 K: x 0.581124 y 0.397741
M2 Planck  2000 K: x 0.524058 y 0.415869
M2 Planck  2856 K: x 0.446241 y 0.408663
M2 Planck  4000 K: x 0.379838 y 0.377370
M2 Planck  6500 K: x 0.313426 y 0.323904
M2 Planck 10000 K: x 0.280726 y 0.288469
M2 Planck 20000 K: x 0.256667 y 0.257776
H4 x strictly decreasing along the locus: True
M2 A formula      : x 0.446275 y 0.408674
M2 skimage A      : x 0.447570 y 0.407450   colour-science A: x 0.447580 y 0.407450
H2 |A formula - colour-science A| 1.305e-03 1.224e-03; - skimage A 1.295e-03 1.224e-03 (bar 0.005)
H3 |Planck 2856 K - A formula| 3.448e-05 1.050e-05 (bar 1e-4)
M3 derived NPM:
 [[0.412391 0.357584 0.180481]
 [0.212639 0.715169 0.072192]
 [0.019331 0.119195 0.950532]]
C4 max|NPM (1,1,1) - D65 XYZ|: 1.110e-16 (tol 1e-12); C4b Y row x1e4: [2126, 7152, 722] (color spaces: [2126, 7152, 722])
H5a max|NPM - skimage xyz_from_rgb|: 3.052e-04 (bar 5e-4)
H5b NPM rounded to 4 decimals equals colour-science matrix: True
M3 whites (X, Y, Z): from xy 0.950456 1.000000 1.089058 | skimage D65 0.950470 1.000000 1.088830 | skimage matrix row sums 0.950456 1.000000 1.088754
H6 |skimage row-sum Z - skimage D65 Z|: 7.600e-05 (bar > 1e-5)
H7 monochromatic stimuli inside the sRGB triangle: 0 of 401; largest min(r,g,b)/(r+g+b): -1.8911e-02
M4 450 nm: x 0.159250 y 0.015480
M4 520 nm: x 0.081028 y 0.820368
M4 580 nm: x 0.513351 y 0.486398
M4 650 nm: x 0.720468 y 0.279532
D1 reflectance min 0.020000 max 0.980000 within [0.02, 0.98]: True
D2 every patch uniform: True ; distinct patches: 24
D3 null-space dimension 5; a_max per row: 0.1800 0.2123 0.7787 0.5872 0.2617 0.4800 all > 0.05: True
L1 6500 K | white x 0.3134 y 0.3239 | max|dXYZ| 1.386e-13 | max dE 5.959e-13 | median dE 0.000000 | identical 8-bit pairs 18 of 18 (t = 1: 6 of 6) | clipped values 0
   row 0 dE at t = 1/3, 2/3, 1: 0.000000 0.000000 0.000000  increasing: False
   row 1 dE at t = 1/3, 2/3, 1: 0.000000 0.000000 0.000000  increasing: False
   row 2 dE at t = 1/3, 2/3, 1: 0.000000 0.000000 0.000000  increasing: False
   row 3 dE at t = 1/3, 2/3, 1: 0.000000 0.000000 0.000000  increasing: False
   row 4 dE at t = 1/3, 2/3, 1: 0.000000 0.000000 0.000000  increasing: False
   row 5 dE at t = 1/3, 2/3, 1: 0.000000 0.000000 0.000000  increasing: False
L2 CIE A | white x 0.4463 y 0.4087 | max|dXYZ| 5.374e+00 | max dE 1.229e+01 | median dE 3.534472 | identical 8-bit pairs 2 of 18 (t = 1: 0 of 6) | clipped values 8192
   row 0 dE at t = 1/3, 2/3, 1: 0.315880 0.632661 0.950351  increasing: True
   row 1 dE at t = 1/3, 2/3, 1: 2.613658 5.264745 7.954711  increasing: True
   row 2 dE at t = 1/3, 2/3, 1: 2.314261 4.653049 7.017614  increasing: True
   row 3 dE at t = 1/3, 2/3, 1: 4.155715 8.250794 12.287674  increasing: True
   row 4 dE at t = 1/3, 2/3, 1: 2.562752 5.075301 7.541533  increasing: True
   row 5 dE at t = 1/3, 2/3, 1: 0.960171 1.931138 2.913229  increasing: True
L3 3-band | white x 0.3306 y 0.3093 | max|dXYZ| 8.742e+00 | max dE 1.973e+01 | median dE 5.786312 | identical 8-bit pairs 0 of 18 (t = 1: 0 of 6) | clipped values 6144
   row 0 dE at t = 1/3, 2/3, 1: 1.628590 3.256170 4.882827  increasing: True
   row 1 dE at t = 1/3, 2/3, 1: 1.409393 2.820184 4.232431  increasing: True
   row 2 dE at t = 1/3, 2/3, 1: 5.234348 10.618848 16.164726  increasing: True
   row 3 dE at t = 1/3, 2/3, 1: 5.354708 11.101455 17.339298  increasing: True
   row 4 dE at t = 1/3, 2/3, 1: 6.388196 12.945188 19.734713  increasing: True
   row 5 dE at t = 1/3, 2/3, 1: 3.106742 6.217915 9.333966  increasing: True
C5 is the L1 max|dXYZ| above (tol 1e-9); C6 max|own Lab - skimage xyz2lab|: 0.000e+00 (tol 1e-10)
H10 median dE: L3 5.786312 > L2 3.534472: True

The checks. All held. The two codings of Wyman et al.‘s fit agree exactly (C1), so Table 1 and Listing 1 describe the same functions; the A formula reproduces colour-science’s docstring to the 7 decimals printed there (C2) and coincides with the ITS-90 Planck shape at 2855.54 K to rounding (C3); the derived matrix maps white to white to rounding (C4) and its luminance row rounds to Color Spaces’ 0.2126, 0.7152 and 0.0722 (C4b); the metamers differ from their bases under L1 by $1.4 \times 10^{-13}$ on a scale where white is 100 (C5); and our CIELAB matches scikit-image’s to the last bit (C6). The data checks held too: the reflectances stay within $[0.02, 0.98]$ and touch both limits, as the amplitude rule intends (D1); every patch is uniform and all 24 are distinct (D2); and the null space has dimension 5, with the smallest amplitude limit 0.18 (D3).

The fit and the white points. H1 held: the fit puts the equal-energy white at $(0.332998, 0.333716)$, within $4 \times 10^{-4}$ of $(1/3, 1/3)$. Since E is exactly $(1/3, 1/3)$ for the tabulated functions, this deviation is error of the fit, together with the range cut at 380 and 780 nm; the output does not separate the two. H2 held: the fit puts A 0.0012 to 0.0013 from the tabulated A in each coordinate, against scikit-image’s value and against colour-science’s alike, which agree with each other to $10^{-5}$. That is the size of error to expect here from the fit and the 380–780 nm range in chromaticity: about $4 \times 10^{-4}$ at E and $1.3 \times 10^{-3}$ at A. H3 held: the ITS-90 blackbody at 2856 K and the CIE formula at 2848 K with the older $c_2$ differ by $3.4 \times 10^{-5}$ in $x$, far below the fit’s error, because their temperatures differ by only 0.46 K once the constants are matched. H4 held: $x$ falls strictly from 1500 K to 20000 K along the locus. The 6500 K blackbody lands at $(0.313426, 0.323904)$, while D65 is $(0.3127, 0.3290)$: 0.0051 apart in $y$ (our arithmetic), about four times the fit’s error at A, so in our reading D65 is not the chromaticity of a 6500 K blackbody, and L1 is a stand-in for daylight, not for D65.

The sRGB matrix. H5a and H5b held. The derived matrix rounds exactly to colour-science’s 4-decimal matrix, which is consistent with that matrix being this derivation rounded to four decimals. scikit-image’s 6-decimal matrix differs from it by up to $3.05 \times 10^{-4}$, so it is not the same derivation at higher precision. H6 held: scikit-image’s matrix sums to a white $Z$ of 1.088754, its D65 constant is 1.08883, and the D65 xy gives 1.089058: two different “D65” whites in one library and a third from the chromaticity, which agree in $Z$ only when rounded to three decimals. Its source says how its constant was chosen: the D65/2° white is copied from an older constant for backward compatibility, while its other whites are computed from chromaticities taken from Wikipedia. In our reading, a test that compares two implementations at $10^{-4}$ has to fix which of these whites it uses.

The spectral locus. H7 held: none of the 401 monochromatic stimuli lies inside the sRGB triangle. The largest value of $\min(r, g, b)/(r + g + b)$ over all wavelengths is $-0.0189$, so even the spectral colour closest to the triangle needs a negative amount of one primary. Read as functions of wavelength, the three rows of NPM⁻¹ applied to $(\bar{x}, \bar{y}, \bar{z})$ are the colour-matching functions of the sRGB primaries (the matching functions of the Formal Definition), and the check says that at every wavelength at least one of them is negative. The 520 nm stimulus, at $(0.081, 0.820)$ under the fit, lies far beyond the sRGB green $(0.300, 0.600)$.

Metamers that match. H8 held: under L1 the largest $\Delta E^*_{ab}$ of the 18 pairs is $6 \times 10^{-13}$, and all 18 pairs get identical 8-bit codes, the first half of H11. The per-row lines under L1 print increasing: False because every $\Delta E^*_{ab}$ there is zero (to six decimals, and $6 \times 10^{-13}$ at most), so no row can increase; H9 is defined for L2 and L3 only.

Metamers that break. The same patches under the other two lights are no longer matched. Under A the median pair differs by $\Delta E^*_{ab} = 3.53$ and the largest by 12.29; under the three-band light the median is 5.79 and the largest 19.73. H9 held: in all six rows and under both lights the colour difference rises strictly with the amplitude of the null-space component, and close to proportionally, as row 3 under A shows (4.16, 8.25, 12.29). The XYZ difference is exactly proportional to the amplitude by construction, and in our reading CIELAB is close enough to linear over these differences for the proportionality to survive. H10 held: the median is larger under the three-band light. It does not hold row by row, which the hypothesis did not claim: row 1, the bump at 480 nm, breaks more under A (7.95) than under the three-band light (4.23), while row 0, the bump at 430 nm, breaks far more under the three-band light (4.88 against 0.95). The output does not explain the difference between rows; in our reading, how much a given pair breaks under a light depends on where its null-space component and its base reflectance lie relative to that light’s spectrum, so neither kind of light breaks every pair more. The second half of H11 held: under both A and the three-band light none of the six largest-amplitude pairs keeps identical 8-bit codes. Two pairs with smaller amplitudes still render identically under A, although every pair under A has $\Delta E^*_{ab} \geq 0.32$; in our reading 8-bit quantization and clipping can hide small differences (the output counts 8192 clipped channel values under A and 6144 under the three-band light, that is 8 and 6 patch channels of 1024 pixels each, against none under L1, since the orange light and the narrow bands push some patches outside the sRGB gamut when no adaptation is applied).

The image therefore makes the textbook statement concrete. An observer that sees only the three numbers under L1 cannot tell any metamer from its base, but a change of light separates them; in this example the three-band light separated them more than A did in the median, though not in every row.

Process note. The plan for this example listed nine departures from the batch plan’s description before the pilot: (1) illuminant A is computed from colour-science’s CIE 015:2004 formula, 2848 K with $c_2 = 1.435 \times 10^{-2}$ m K, and printed beside the ITS-90 Planck spectrum at 2856 K; the planning research had proposed evaluating Planck at “2856 × 1.4388/1.435 ≈ 2863.6 K”, which is wrong, since the formula’s own temperature is 2848 K and the ITS-90 equivalent is 2855.54 K (C3); (2) A is also compared with colour-science’s tabulated value, not only scikit-image’s; (3) the matrix is also compared with colour-science’s 4-decimal matrix (H5b) and its luminance row with Color Spaces’ weights (C4b); (4) the metamers form an amplitude series under two breaking lights, the smooth A and a three-band light (H9, H10), where the batch plan named one light; (5) the matching light is a 6500 K blackbody, because D65’s tabulated spectrum is data that is not downloaded; (6) $\Delta E^*_{ab}$ with our CIELAB cross-checked against scikit-image (C6) and the identical-code count (H11) were added as measures; (7) checks C1–C6 and data checks D1–D3 were added; (8) a gamut fraction of random smooth reflectances proposed at planning was dropped, because its generator repeats one in the color constancy example; (9) data check D2 was repaired before the pilot. One planning script computed the checks’ reachable values: it ran once, after the hypotheses and the checks’ bars were written, and printed only the checks’ values. In it, D2’s first form, which tested each patch’s uniformity as “patch std == 0”, failed: it printed False while the distinct-patch count was 24, because the standard deviation of identical floating-point values is not exactly zero. In the plan’s words, “D2’s uniformity test is changed before the pilot to “every pixel of a patch equals its first pixel exactly” (np.all(p == p[0]))“. The bars of C1–C6, D2, D3, H1–H7 and H9–H11 were written before any script ran, and these checks and hypotheses are tests. Two are illustrations: D1, whose slack of $10^{-12}$ around $[0.02, 0.98]$ was added after the planning script printed the reflectance range as 0.020000 to 0.980000; and H8, written before any script, whose outcome the planning script had already shown in its XYZ form (C5) before the pilot. The pilot’s output was identical to the first run’s, and no setting, measure or printed value changed after the pilot. After the first run, one code comment was reworded (a name for the null-space component was replaced by a description); the printed output was unchanged.

Common Misconceptions

• “XYZ values are amounts of three real lights.” They are coordinates in a space whose basis lights do not exist; the CIE functions were chosen to be non-negative at the cost of having no physical primaries (visionbook, §8.3.6).
• “The chromaticity diagram is a map of perceived colour differences.” It is a projection in which mixtures lie on lines. Equal distances in it are not equal visible differences, and its gamut areas are not perceptual measures; CIELAB and its colour differences exist for that purpose (color spaces).
• “Equal XYZ means equal colour under any conditions.” It means a match for the standard observer under the light in the calculation. The example’s pairs match under one light to $10^{-13}$ and differ by up to 12 and 20 $\Delta E^*_{ab}$ under two others.
• “D65 is the colour of a 6500 K blackbody.” In the example the 6500 K blackbody lies 0.005 from D65 in $y$, more than the fit’s error; D65 belongs to the CIE daylight series, whose spectra colour-science’s sd_CIE_illuminant_D_series computes from tabulated basis functions, not from Planck’s law; that function’s docstring example takes D65’s nominal temperature as $6500 \times 1.4388/1.4380$ K, in our reading another correction for a changed $c_2$.
• “A white point is a property of the light alone.” It also depends on the observer: D65 is $(0.31270, 0.32900)$ for the 1931 observer and $(0.31382, 0.33100)$ for the 1964 observer.
• “Every library’s sRGB matrix and D65 are the same numbers.” They agree to about three or four decimals. colour-science’s 4-decimal matrix equals the derivation rounded to four decimals; scikit-image’s 6-decimal matrix differs from the derivation by $3 \times 10^{-4}$ and does not reproduce its own D65 constant.
• “Illuminant A is a 2856 K blackbody with today’s constants.” Its formula uses 2848 K with the older $c_2 = 1.435 \times 10^{-2}$ m K, equivalent to 2855.54 K with the ITS-90 constant; the difference from 2856 K is negligible in chromaticity, but it matters in a test that compares spectra.
• “A camera sees colorimetric values.” Its filters do not match the CIE functions (image formation), so in our reading it has its own metamers; the colour correction matrix that maps camera RGB to a standard space (digital camera) cannot separate two surfaces that the camera’s RGB merged, nor, being invertible, merge two that it separated.

Where It Is Used

Color spaces defines sRGB, CIELAB and related spaces through CIE XYZ, primaries and white points, and gamma correction applies their transfer functions and CIE lightness. In a camera, the image sensor and its Bayer filter sample the spectrum through filters that are not the CIE functions, and the digital camera pipeline maps the result to a colorimetric space with white-balance gains and a colour correction matrix; image formation follows the spectral chain from light to pixel. Color constancy estimates the light whose white point colorimetry assumes known, and radiometry and photometry shares the luminous weighting of $Y$ and the Planck spectra used for illuminants. In our reading, metamerism is also why a spectral or multispectral measurement can distinguish materials that an RGB image, or a person, cannot: more bands leave a smaller null space.