PL. 127 · COLOUR / COLOUR MATCHING / NULL SPACE
Metameric Blacks
Günter Wyszecki, 1953 · the colour matching fit: Wyman, Sloan & Shirley, 2013
OPEN THE LIVE PLATE ▸DEFINITION
A = [x̄(λ_k) ; ȳ(λ_k) ; z̄(λ_k)], a 3 × N matrix b ∈ null A, dim null A = N − 3 the metameric blacks s± = s ± α b ⇒ A s+ = A s− , α at the bound s ± α b ≥ 0
NOTES
Every other plate in this order draws a difference. This one draws an equality. A spectrum arriving at the eye is a long list of numbers, one per wavelength; what leaves is three, because there are three kinds of cone and no more. Sample the spectrum at N wavelengths and that reduction is a matrix — three rows, N columns, each row one of the colour matching functions. Any matrix that is wider than it is tall has a null space, and this one has dimension N minus three: an enormous room full of spectral shapes the eye maps to exactly zero. Wyszecki named them metameric blacks in 1953. Add one to any spectrum and the colour does not move. The two curves above are that addition and that subtraction, one base spectrum plus and minus the same black, and the band below is the colour, integrated twice over — the left half from one of the two and the right half from the other, each integrated over its own curve and neither copying the other. There is a join down the middle of it and no way to find it. Turn the amplitude up and the curves pull apart until one of them touches the axis: that is the only wall in the room, because a spectrum cannot carry negative light, and where it touches is where the base spectrum was dimmest. Turn the sample count down and the room shrinks, since the null space is N minus three dimensional and three wavelengths leave nothing at all. The drift is a rotation inside that space — the black turns through a plane of itself, the pair writhes, and the colour never moves, because it cannot.
PROVENANCE
- Origin
- G. Wyszecki, "Valenzmetrische Untersuchung des Zusammenhanges zwischen normaler und anomaler Trichromasie", Die Farbe 2, 1953, 39–52 — the paper in which a radiometric function is first decomposed into a fundamental part and a residue the eye cannot see, and the reason that residue is called a Wyszecki black. Die Farbe was never retro-registered and no DOI resolves to it, so this is a citation rather than a link; the volume and the opening page are the ones the later literature agrees on, Cohen and Kappauf below among them.
- The later formalisation
- J. B. Cohen and W. E. Kappauf, "Metameric Color Stimuli, Fundamental Metamers, and Wyszeckis Metameric Blacks", The American Journal of Psychology 95(4), 1982, 537–564 — the src link above, and the paper that turned the idea into linear algebra. Their matrix R = A(AᵀA)⁻¹Aᵀ is the orthogonal projector onto the row space of the colour matching matrix, so I − R projects onto exactly the blacks. This plate builds the same projector by Gram–Schmidt instead of by inversion, which is the same operator arrived at more stably. Their fuller treatment, "Color Mixture and Fundamental Metamers: Theory, Algebra, Geometry, Application", is a separate and later paper, Am. J. Psychol. 98(2), 1985, 171–259; the two are easy to conflate and are not the same work.
- Why the map is a matrix at all
- Hermann Grassmann, "Zur Theorie der Farbenmischung", Annalen der Physik 165(5), 1853, 69–84 (doi:10.1002/andp.18531650505). Nothing here is possible unless colour matching is linear — unless mixing two lights adds their matches. Grassmann is where that was set out as law rather than assumed, and it is the whole reason a null space exists to draw.
- The colour matching functions, and why the tables are absent
- Turning a spectrum into three numbers needs the CIE 1931 standard observer, whose tabulated x̄ ȳ z̄ values are asserted-copyright and are therefore not in this tree and never will be. What is here instead is the published analytic fit: C. Wyman, P.-P. Sloan and P. Shirley, "Simple Analytic Approximations to the CIE XYZ Color Matching Functions", Journal of Computer Graphics Techniques 2(2), 2013, 1–11 (jcgt.org/published/0002/02/01/) — seven piecewise-continuous Gaussian lobes, three for x̄ and two each for ȳ and z̄, evaluated here from the amplitudes, centres and reciprocal widths of Equation 4 and Table 1. Those coefficients are the content the paper exists to publish, offered so that people implement them instead of copying tables; using them is attribution, given here in full, and not a licence. No code was taken: the paper ships a C listing and a supplemental repository and neither was read into this implementation. Blackbody, dispersion and gamut carry the same fit, written out separately in each — four copies of seven lines rather than one shared helper, because a specimen has to be readable on its own in the drawer.
- Standing
- Public domain. The construction is Wyszeckis decomposition and Cohen and Kappaufs projector, implemented from the linear algebra they set out; the observer is the analytic fit named above. Nothing is taken from any colour science or rendering codebase.
- The null space, constructed rather than solved
- No numerical library and no matrix inversion. The three rows of the matrix are Gram–Schmidt orthonormalised against each other in closed form, working in coefficient space over (x̄, ȳ, z̄) so that the inner product is one 3 × 3 form built in the same pass that integrates the spectrum. Two test functions — cosines of nulldir and nulldir + 1 half cycles across the band — are then projected off that basis, and what is left is orthogonal to all three rows, which is exactly what it means to be in the null space. The two leftovers are orthonormalised against each other in turn, so the drift is a genuine rotation of a unit circle inside the space rather than a wander through an oblique pair. Everything reduces to five scalars: two coefficients on the test functions and a triple on the matching functions themselves, which is why the shader can evaluate the black at any wavelength without an array.
- Where the amplitude stops
- A metameric black is signed, so adding a large enough multiple of one drives the spectrum negative, and a spectrum with negative light in it is not a spectrum. That is the real and only bound on this plate, and the amplitude dial names a position inside it rather than an absolute size — 1 is exactly where the lower curve first touches zero, which is dispersions transmitting-window trick applied to a different constraint. It is also why every value of the dial renders: the bound is recomputed each frame from the spectrum actually on screen, so there is no setting at which the pair silently leaves physics. The wall moves as the black rotates, and the wavelength at which the curve kisses the axis moves with it.
- What the equality is relative to
- Exactly this observer, sampled at exactly these N wavelengths. That is not a hedge, it is what metamerism is: two spectra match for one eye and one sampling, and change either and the match can break. The plate is honest about this in the one way that matters — the colour of each half of the band is integrated over the same N samples the black was built from, so the two agree because the mathematics says so, not because the plate drew one rectangle twice. Move the sample count and both the black and the colour change together. The curve drawn between the sample dots is a straight-line interpolation and is scenery: the object is the dots, and only the dots are constrained.
- Checked, not just plotted
- Three independent figures. First, the equality: taken through a float32 emulation of the shaders own arithmetic across the whole declared box — 101,824 constant tuples — the two halves of the band differ by at most 0.00055 of one level in 255 before rounding, and 6 of those 101,824 land close enough to a rounding boundary that Math.round sends them to adjacent levels. Everywhere else the two halves are byte-identical, which is the strongest form the claim can take on an 8-bit screen. Second, the observer: an equal-energy spectrum is CIE illuminant E, whose chromaticity is x = y = 1/3 by definition — a pair of numbers nothing in this implementation has seen. Setting the tilt to zero and integrating gives x = 0.33304, y = 0.33376 at 36 samples, 0.00052 away in the chromaticity plane, and that residual is the price of the fit plus the cut of the band at 780 nm rather than 830; it is already 0.00051 at 64 samples and does not move again out to 4096. Third, the gamut: every colour this plate can mix was checked against sRGB across the whole tilt and sample box and none of them is outside it, so unlike blackbody there is nothing here being clamped, and the band is the colour rather than the nearest one a screen has.
- Constants
- amp is the amplitude as a fraction of the non-negativity bound, described above. tilt shapes the base spectrum: it is exp(tilt·(u − ½)) with u running 0 to 1 across the band and the peak normalised to 1, so the family is the simplest smooth positive one-parameter family there is, running from a cool blue-white at −1.8 through equal energy at 0 to a warm amber at +1.8. Its range is where it is because the bound above is set by wherever the base spectrum is dimmest, and a steeper tilt is a dimmer far end with less room for a black in it. Measured, the narrowest the two curves ever come at full amplitude — worst over all eight null directions and every angle of the drift — is 0.34 of the base peak at ±1.8, and half of that, 0.17, at ±2.5, which is why the range stops where it does. nulldir picks the two test functions and so which plane of the null space the drift turns in — low values are slow rolling shapes, high ones are fine ripples. Its ceiling and the sample floor are one decision: cos of nulldir + 1 half cycles has to be resolvable by N samples, and at N equal to the half-cycle count the test function is zero at every sample and the plane collapses. Measured, that happens at 6 and 8 samples for the upper half of this dial and nowhere at or above 12, so the sample floor is 12 and the smallest ratio between the two null directions across the whole declared box is 0.81 — a plane, never a line. samples is locked from regenerate for the reasons in taxonomy.js. Liveness: 319 of 320 random four-constant tuples drawn uniformly from the whole declared box separate the two curves by more than a tenth of the base spectrum peak at every angle of the drift, and the single miss reaches 0.094 rather than 0 — two curves about ten pixels apart on a card, not a blank plate. All 40 tuples drawn the way regenerate actually jitters are live with a worst separation of 0.399, and no slider has a dead position when travelled alone.
- Source
- doi:10.2307/1422186
TOUCHDESIGNER · GLSL
The same shader this plate runs, reframed for a GLSL TOP. Pasted bare it renders the published constants as a still frame; wire absTime.seconds into u_t on the Vectors page to animate it.
// FORMA — PL. 127 · METAMERIC BLACKS — Günter Wyszecki, 1953 · the colour matching fit: Wyman, Sloan & Shirley, 2013
// A = [x̄(λ_k) ; ȳ(λ_k) ; z̄(λ_k)], a 3 × N matrix
// b ∈ null A, dim null A = N − 3 the metameric blacks
// s± = s ± α b ⇒ A s+ = A s− , α at the bound s ± α b ≥ 0
// TouchDesigner port — paste into a GLSL TOP's pixel shader. Set the
// resolution on the TOP's Common page. As pasted it renders the published
// constants as a still frame; to animate, add a uniform named u_t on the
// GLSL TOP's Vectors 1 page with the expression absTime.seconds.
// Constants are consts — edit to tweak; comments give the measured range.
// Written from the published mathematics, not adapted from any code.
#define u_res (uTDOutputInfo.res.zw)
uniform float u_t; // absTime.seconds on the Vectors page; unset = still
const float u_phase = 0.9622; // this plate's own grid phase, 0..1
// FORMA's COLOUR accent as cosine-gradient coefficients
const vec3 u_pal_a = vec3(0.46, 0.2002, 0.3896);
const vec3 u_pal_b = vec3(0.5, 0.2176, 0.4235);
const vec3 u_pal_c = vec3(1, 1, 1);
const vec3 u_pal_d = vec3(0, 0.05, 0.1);
const float p_amp = 0.8; // black amplitude (1 = the non-negativity bound) · live 0.25 .. 1
const float p_tilt = 0.9; // base spectrum tilt · live -1.8 .. 1.8
const float p_nulldir = 3.0; // null direction (half cycles across the band) · live 1 .. 8
const float p_samples = 36.0; // spectral samples · live 12 .. 88
/* The order's ramp — the same cosine formulation the JS kit uses, so a
plate keeps its classification colour in either language. */
vec3 ramp(float t){
return clamp(u_pal_a + u_pal_b * cos(6.28318530718 * (u_pal_c * t + u_pal_d)), 0.0, 1.0);
}
/* Sawtooth and triangle on this plate's phase, mirroring the JS kit. */
float cycle(float t, float period){ return fract(t / period + u_phase); }
float pingpong(float t, float period){
float u = cycle(t, period);
return u < 0.5 ? u * 2.0 : 2.0 - u * 2.0;
}
/* The colour matching map is linear (Grassmann, 1853), so sampling a spectrum
at N wavelengths makes it a 3 x N matrix, and a matrix wider than it is tall
has a null space. Its elements are Wyszecki's metameric blacks: spectral
shapes the eye maps to exactly zero. This shader builds the matrix, the null
space and the two spectra per fragment, because a fragment shader has
nowhere to put a result — the numbers are identical to the JS path, so the
two are one picture sampled twice.
The CIE tabulations are asserted-copyright and are not in this tree; what is
here is the published analytic fit of Wyman, Sloan and Shirley (JCGT 2(2),
2013), evaluated from the coefficients of its own Equation 4 and Table 1.
Every numeric constant below is declared inside the function that uses it,
never at file scope. A composition emits two copies of a body and namespaces
only the function declarations, so a file-scope const would be declared
twice by metamer x metamer and fail to link. */
/* One piecewise-continuous Gaussian lobe: the width differs either side of the
centre, which is how the fit follows an asymmetry a plain Gaussian cannot.
g and d are reciprocal widths. */
float meta_lobe(float lambda, float a, float b, float g, float d){
float x = (lambda - b) * (lambda < b ? g : d);
return a * exp(-0.5 * x * x);
}
/* x-bar, y-bar, z-bar of the CIE 1931 two-degree observer: three lobes, then
two and two. The negative amplitude is the real dip in x-bar near 500 nm. */
vec3 meta_cmf(float lambda){
return vec3(
meta_lobe(lambda, 0.362, 442.0, 0.0624, 0.0374)
+ meta_lobe(lambda, 1.056, 599.8, 0.0264, 0.0323)
+ meta_lobe(lambda, -0.065, 501.1, 0.0490, 0.0382),
meta_lobe(lambda, 0.821, 568.8, 0.0213, 0.0247)
+ meta_lobe(lambda, 0.286, 530.9, 0.0613, 0.0322),
meta_lobe(lambda, 1.217, 437.0, 0.0845, 0.0278)
+ meta_lobe(lambda, 0.681, 459.0, 0.0385, 0.0725));
}
/* CIE XYZ to linear sRGB (IEC 61966-2-1 primaries at D65), then that
standard's transfer function — a linear toe under a shifted power curve. */
vec3 meta_rgb(vec3 c){
return vec3( 3.2404542 * c.x - 1.5371385 * c.y - 0.4985314 * c.z,
-0.9692660 * c.x + 1.8760108 * c.y + 0.0415560 * c.z,
0.0556434 * c.x - 0.2040259 * c.y + 1.0572252 * c.z);
}
vec3 meta_encode(vec3 c){
c = clamp(c, 0.0, 1.0);
return mix(12.92 * c, 1.055 * pow(c, vec3(1.0 / 2.4)) - 0.055, step(0.0031308, c));
}
/* An exponential tilt across the band, peak normalised to 1. Strictly
positive everywhere: the amplitude bound is set by wherever this is
smallest, and a base that touched zero would admit no black at all. */
float meta_base(float u, float tilt){
return exp(tilt * (u - 0.5) - abs(tilt) * 0.5);
}
/* The discrete inner product of two vectors written in coefficient space over
(x-bar, y-bar, z-bar): everything the orthogonalisation needs to know about
the observer is the 3 x 3 Gram matrix of those three rows. */
float meta_ip(vec3 a, vec3 b, mat3 G){ return dot(a, G * b); }
/* The metameric black at one wavelength. Five scalars are its whole
description — two coefficients on the test functions and a triple on the
matching functions — which is what lets a fragment shader evaluate it
anywhere without an array. The matching functions arrive already evaluated,
because every caller here has just needed them for something else and they
are seven exponentials each. */
float meta_black(float u, vec3 c, float dir, float A1, float A2, vec3 Q){
const float PI = 3.141592653589793;
return A1 * cos(dir * PI * u) + A2 * cos((dir + 1.0) * PI * u) + dot(Q, c);
}
/* Chromaticity, not radiance: the base level of a spectrum is arbitrary here,
so each half of the band is scaled to the brightest sRGB triple carrying its
own colour. Measured across the whole declared box, nothing this plate can
mix is outside the gamut, so the max() is a formality rather than a clamp. */
vec3 meta_swatch(vec3 xyz){
vec3 lin = meta_rgb(xyz / max(xyz.y, 1e-20));
return meta_encode(max(lin, 0.0) / max(max(lin.r, lin.g), max(lin.b, 1e-9)));
}
/* 1 inside a mark and 0 outside it, with the transition spread over the
outer part of its width so the edge is antialiased rather than stepped. */
float meta_mark(float e0, float e1, float x){ return 1.0 - smoothstep(e0, e1, x); }
vec3 plate(vec2 uv){
/* The plate in three registers, as fractions of its height. */
const float PLOT_T = 0.055, PLOT_B = 0.648;
const float RULE_T = 0.674, RULE_B = 0.708;
const float BAND_T = 0.756, BAND_B = 0.960;
const float TICK_A = 0.722, TICK_B = 0.744, TICK_C = 0.972, TICK_D = 0.994;
const float LO = 380.0, SPAN = 400.0;
/* The largest sRGB channel the matching functions reach when taken straight
through the matrix — measured over the band, at 605.7 nm. */
const float HUE_PEAK = 2.503;
const float PI = 3.141592653589793;
float asp = u_res.y / u_res.x;
/* Four is the floor the mathematics has: a 3 x N matrix has no null space
at all until N exceeds 3. The slider starts at 12 for a separate reason
(see the header), so this max is never the binding one — it is here so
that neither path can be handed a count that makes the object vanish. */
int n = max(4, int(p_samples + 0.5));
float fn = float(n);
float dir = floor(p_nulldir + 0.5);
/* ---- pass one: the observer's Gram matrix, the test functions against it,
and the base spectrum's own tristimulus. Constant loop bound at the
slider's own ceiling, break on the live value, mandelbrot fashion. */
mat3 G = mat3(0.0);
vec3 d1 = vec3(0.0), d2 = vec3(0.0), Sxyz = vec3(0.0), ff = vec3(0.0);
for (int k = 0; k < 88; k++){
if (k >= n) break;
float u = (float(k) + 0.5) / fn;
vec3 c = meta_cmf(LO + u * SPAN);
float s = meta_base(u, p_tilt);
float g1 = cos(dir * PI * u), g2 = cos((dir + 1.0) * PI * u);
G[0] += c.x * c; G[1] += c.y * c; G[2] += c.z * c;
d1 += g1 * c; d2 += g2 * c;
ff += vec3(g1 * g1, g1 * g2, g2 * g2);
Sxyz += s * c;
}
/* Gram–Schmidt on the three rows, in coefficient space rather than over the
N samples: a vector is three numbers and the inner product is G. Cohen and
Kappauf write the same projector as A(A'A)inv A'; this is that operator
built without an inversion. */
vec3 e1 = vec3(1.0, 0.0, 0.0);
e1 = e1 / sqrt(meta_ip(e1, e1, G));
vec3 v2 = vec3(0.0, 1.0, 0.0);
v2 = v2 - e1 * meta_ip(v2, e1, G);
vec3 e2 = v2 / sqrt(meta_ip(v2, v2, G));
vec3 v3 = vec3(0.0, 0.0, 1.0);
v3 = v3 - e1 * meta_ip(v3, e1, G);
v3 = v3 - e2 * meta_ip(v3, e2, G);
vec3 e3 = v3 / sqrt(meta_ip(v3, v3, G));
/* Projecting the two test functions off that basis leaves two vectors
orthogonal to all three rows, which is what being in the null space means.
Their mutual inner products come from numbers already in hand. */
vec3 a1 = vec3(dot(d1, e1), dot(d1, e2), dot(d1, e3));
vec3 a2 = vec3(dot(d2, e1), dot(d2, e2), dot(d2, e3));
vec3 q1 = a1.x * e1 + a1.y * e2 + a1.z * e3;
vec3 q2 = a2.x * e1 + a2.y * e2 + a2.z * e3;
float b11 = ff.x - dot(a1, a1);
float b12 = ff.y - dot(a1, a2);
float b22 = ff.z - dot(a2, a2);
float n1 = sqrt(max(b11, 1e-12));
/* Orthonormalising the second against the first makes the drift a rotation
of a circle rather than a wander round an oblique pair. The floor guards
the case where the sample count cannot resolve the two test functions
apart, which is why the sample floor is 12; measured, the smallest this
ratio reaches inside the declared box is 0.81. */
float n2 = max(sqrt(max(b22 - b12 * b12 / (n1 * n1), 0.0)), 1e-4 * n1);
float mu = b12 / (n1 * n1);
/* The drift: a rotation inside the null space. cycle() carries this plate's
own PHASE, so no two cards sit at the same angle. */
float th = 6.28318530718 * cycle(u_t, 40.0);
float A1 = cos(th) / n1 - sin(th) * mu / n2;
float A2 = sin(th) / n2;
vec3 Q = -(A1 * q1 + A2 * q2);
/* ---- pass two: the bound, and the black's own tristimulus. That second
sum is the claim, computed rather than assumed — if the construction is
right it is zero, and the two halves of the band are built by adding and
subtracting it. */
float limit = 1e20;
vec3 D = vec3(0.0);
for (int k = 0; k < 88; k++){
if (k >= n) break;
float u = (float(k) + 0.5) / fn;
vec3 c = meta_cmf(LO + u * SPAN);
float b = meta_black(u, c, dir, A1, A2, Q);
/* A spectrum cannot carry negative light, so the multiple that may be
added is capped by whichever sample runs out of base spectrum first. */
if (abs(b) > 1e-9) limit = min(limit, meta_base(u, p_tilt) / abs(b));
D += b * c;
}
float scale = p_amp * limit;
float PMAX = 1.0 + p_amp; // the base peak is 1, so this is the ceiling
/* ---- the exhibit ---------------------------------------------------- */
if (uv.y > BAND_T && uv.y < BAND_B){
/* Each half from its own spectrum. Factoring the constant out of the sum
is arithmetic, not a shortcut: sum of (s +- scale.b).cbar really is
sum of s.cbar +- scale times sum of b.cbar, and that second sum is the
quantity just measured rather than one assumed. The join is at 0.5 and
nothing marks it. */
return meta_swatch(uv.x < 0.5 ? Sxyz + scale * D : Sxyz - scale * D);
}
if (uv.y > RULE_T && uv.y < RULE_B){
/* The wavelength rule, so the horizontal reads as wavelength: the matching
functions straight through the sRGB matrix — the fit used only for what
it was fitted to, the absolute values — with negatives clipped, because
a single wavelength is outside the gamut nearly everywhere. */
vec3 l = meta_rgb(meta_cmf(LO + uv.x * SPAN));
return meta_encode(max(l, 0.0) / HUE_PEAK);
}
if ((uv.y > TICK_A && uv.y < TICK_B) || (uv.y > TICK_C && uv.y < TICK_D)){
/* Two short ticks, above the band and below it, aimed at the join. */
return vec3(0.016, 0.024, 0.039)
+ vec3(0.55, 0.58, 0.66) * meta_mark(0.0018, 0.0034, abs(uv.x - 0.5));
}
if (uv.y > PLOT_T - 0.02 && uv.y < RULE_T){
vec3 col = vec3(4.0, 6.0, 10.0) / 255.0;
/* The zero axis, so that a curve touching it is legible as touching. */
col = mix(col, vec3(0.28, 0.30, 0.36), meta_mark(0.0030, 0.0052, abs(uv.y - PLOT_B) * asp));
/* The two spectra between the outer samples, straight-line interpolated
between them. The object is the dots; the line between them is scenery,
because only the samples are constrained. */
float x = uv.x * fn - 0.5;
if (x >= 0.0 && x <= fn - 1.0){
float k0f = clamp(floor(x), 0.0, fn - 2.0);
float fr = clamp(x - k0f, 0.0, 1.0);
float u0 = (k0f + 0.5) / fn, u1 = (k0f + 1.5) / fn;
float s0 = meta_base(u0, p_tilt), s1 = meta_base(u1, p_tilt);
float b0 = meta_black(u0, meta_cmf(LO + u0 * SPAN), dir, A1, A2, Q);
float b1 = meta_black(u1, meta_cmf(LO + u1 * SPAN), dir, A1, A2, Q);
/* The three curves at the two ends of this segment, packed so that one
line of arithmetic serves all of them: x is the base plus the black,
y the base minus it, z the base itself. */
vec3 y0 = vec3(s0 + scale * b0, s0 - scale * b0, s0);
vec3 y1 = vec3(s1 + scale * b1, s1 - scale * b1, s1);
vec3 c0 = PLOT_B - (y0 / PMAX) * (PLOT_B - PLOT_T);
vec3 c1 = PLOT_B - (y1 / PMAX) * (PLOT_B - PLOT_T);
/* A line drawn by vertical distance alone thins wherever it climbs, so
the vertical offset is turned into a perpendicular one through the
segment's own screen slope. */
vec3 slope = (c1 - c0) * fn * asp;
vec3 per = inversesqrt(1.0 + slope * slope);
vec3 vc = mix(c0, c1, fr);
/* Sample dots, on the nearest sample's own centre — the object being
drawn is those points, and the line between them is scenery. */
float kd = clamp(floor(uv.x * fn), 0.0, fn - 1.0);
float ud = (kd + 0.5) / fn;
float sd = meta_base(ud, p_tilt);
float bd = meta_black(ud, meta_cmf(LO + ud * SPAN), dir, A1, A2, Q);
vec2 dv = PLOT_B - (vec2(sd + scale * bd, sd - scale * bd) / PMAX) * (PLOT_B - PLOT_T);
float du = uv.x - ud;
float LW = 0.0075; // line half width, in frame widths
float DR = min(0.0135, 0.45 / fn); // dot radius, same units
col = mix(col, vec3(0.34, 0.36, 0.42),
meta_mark(LW * 0.55, LW, abs(uv.y - vc.z) * asp * per.z) * 0.55);
col = mix(col, vec3(0.80, 0.85, 0.94),
max(meta_mark(LW * 0.55, LW, abs(uv.y - vc.y) * asp * per.y),
meta_mark(DR * 0.6, DR, length(vec2(du, (uv.y - dv.y) * asp)))));
/* The order accent, in the ramp's bright lobe: the ramp is black around
0.5 and these are sparse strokes on a dark ground. */
col = mix(col, ramp(0.92),
max(meta_mark(LW * 0.55, LW, abs(uv.y - vc.x) * asp * per.x),
meta_mark(DR * 0.6, DR, length(vec2(du, (uv.y - dv.x) * asp)))));
}
return col;
}
return vec3(4.0, 6.0, 10.0) / 255.0;
}
out vec4 fragColor;
void main(){
// FORMA's uv runs y-down, matching its canvas; TD's vUV runs up
vec2 uv = vec2(vUV.s, 1.0 - vUV.t);
fragColor = TDOutputSwizzle(vec4(plate(uv), 1.0));
}
NUKE · BLINKSCRIPT
The same shader this plate runs, transpiled to a BlinkScript kernel. Paste it into a BlinkScript node's Kernel Source and press Recompile; every constant arrives as a knob at its published value, and u_t animates with the expression frame/24. Compiled and rendered in Nuke 17.1, then compared against this plate on the page.
// FORMA — PL. 127 · METAMERIC BLACKS — Günter Wyszecki, 1953 · the colour matching fit: Wyman, Sloan & Shirley, 2013
// A = [x̄(λ_k) ; ȳ(λ_k) ; z̄(λ_k)], a 3 × N matrix
// b ∈ null A, dim null A = N − 3 the metameric blacks
// s± = s ± α b ⇒ A s+ = A s− , α at the bound s ± α b ≥ 0
// Nuke port — a BlinkScript kernel. Paste into a BlinkScript node's Kernel
// Source and press Recompile. Every constant arrives as a knob at its published
// value (the comment gives the measured range); u_t is a knob too — animate it
// with the expression frame/24 or leave it at 0 for the still frame. Written
// from the published mathematics, not adapted from any code.
// Transpiled from the shader this plate runs on the page (GLSL ES 3.00):
// vec → float2/3/4, swizzles expanded, GLSL builtins Blink lacks written out
// as forma_ functions, float literals suffixed. Compiled and rendered in a
// real Nuke (17.1v1) and compared against this plate on the page: 34 of 34.
//
// plate() and its helpers are written to a single exit — the loop that runs
// once. That is not a style: Blink 17.1 drops a conditional early return from
// a called function while Vectorize is on, which is the node default, with no
// warning and no error. Written this way it paints correctly as pasted.
kernel Forma_metamer : ImageComputationKernel<ePixelWise>
{
Image<eWrite> dst;
param:
float u_t; // seconds; 0 is the still frame
float p_amp; // black amplitude (1 = the non-negativity bound) · live 0.25 .. 1
float p_tilt; // base spectrum tilt · live -1.8 .. 1.8
float p_nulldir; // null direction (half cycles across the band) · live 1 .. 8
float p_samples; // spectral samples · live 12 .. 88
local:
float2 u_res;
float u_phase;
float3 u_pal_a, u_pal_b, u_pal_c, u_pal_d;
void define(){
defineParam(u_t, "u_t", 0.0f);
defineParam(p_amp, "p_amp", 0.8f);
defineParam(p_tilt, "p_tilt", 0.9f);
defineParam(p_nulldir, "p_nulldir", 3.0f);
defineParam(p_samples, "p_samples", 36.0f);
}
void init(){
u_res = float2(float(dst.bounds.width()), float(dst.bounds.height()));
u_phase = 0.9622f; // this plate's own grid phase, 0..1
// FORMA's COLOUR accent as cosine-gradient coefficients
u_pal_a = float3(0.46f, 0.2002f, 0.3896f);
u_pal_b = float3(0.5f, 0.2176f, 0.4235f);
u_pal_c = float3(1.0f, 1.0f, 1.0f);
u_pal_d = float3(0.0f, 0.05f, 0.1f);
}
/* GLSL builtins Blink lacks, written as templates rather than overload sets.
Blink's operators return expression templates (Swizzle<float,N>), so a call
passing an expression cannot resolve against an overload set on float2
against float3 — measured in Nuke 17.1: a float2 expression is ambiguous
between the two, while scalar-against-vector resolves. A template deduces
the expression's own type, so the ambiguity cannot arise. */
template <class T> T forma_fract(T v){ return v - floor(v); }
template <class T, class S> T forma_mod(T x, S y){ return x - y * floor(x / y); }
/* Blink's own min/max/clamp take no scalar bound against a vector, which GLSL
does; v * 0.0f + b is that bound at the vector's own width, and collapses to
b when v is a scalar, so one template serves both. */
template <class T, class S> T forma_min(T a, S b){ return min(a, a * 0.0f + b); }
template <class T, class S> T forma_max(T a, S b){ return max(a, a * 0.0f + b); }
template <class T, class S> T forma_clamp(T v, S lo, S hi){ return clamp(v, v * 0.0f + lo, v * 0.0f + hi); }
int forma_min(int a, int b){ return min(a, b); }
int forma_max(int a, int b){ return max(a, b); }
/* GLSL step(edge, x) is 1 where x >= edge; floor(sign(x - e) * 0.5 + 1) is
that exactly, equality included, out of builtins Blink does have. */
template <class T, class S> T forma_step(S e, T x){ return floor(sign(x - e) * 0.5f + 1.0f); }
template <class T, class S> T forma_smoothstep(S a, S b, T x){
T t = forma_clamp((x - a) / (b - a), 0.0f, 1.0f);
return t * t * (3.0f - 2.0f * t);
}
template <class T> float forma_distance(T a, T b){ return length(a - b); }
float forma_tanh(float x){ float e = exp(2.0f * x); return (e - 1.0f) / (e + 1.0f); }
float forma_radians(float d){ return d * 0.01745329252f; }
// the page's hash2 is exact uint32; Blink has int, so the shifts are made
// logical by masking and the read-back is lifted into 0 .. 2^32
/* A uint read back as a float. Blink has no unsigned type, so a value past
2^31 arrives as a negative int and float() of it is negative. Measured on
gabor, whose own generator then returned uniforms in [-0.5, 0.5) and drew
a different picture — it compiled, it rendered, and only comparing it with
/* The order's ramp — the same cosine formulation the JS kit uses, so a
plate keeps its classification colour in either language. */
float3 ramp(float t){
return forma_clamp(u_pal_a + u_pal_b * cos(6.28318530718f * (u_pal_c * t + u_pal_d)), 0.0f, 1.0f);
}
/* Sawtooth and triangle on this plate's phase, mirroring the JS kit. */
float cycle(float t, float period){ return forma_fract(t / period + u_phase); }
float pingpong(float t, float period){
float u = cycle(t, period);
return u < 0.5f ? u * 2.0f : 2.0f - u * 2.0f;
}
/* The colour matching map is linear (Grassmann, 1853), so sampling a spectrum
at N wavelengths makes it a 3 x N matrix, and a matrix wider than it is tall
has a null space. Its elements are Wyszecki's metameric blacks: spectral
shapes the eye maps to exactly zero. This shader builds the matrix, the null
space and the two spectra per fragment, because a fragment shader has
nowhere to put a result — the numbers are identical to the JS path, so the
two are one picture sampled twice.
The CIE tabulations are asserted-copyright and are not in this tree; what is
here is the published analytic fit of Wyman, Sloan and Shirley (JCGT 2(2),
2013), evaluated from the coefficients of its own Equation 4 and Table 1.f
Every numeric constant below is declared inside the function that uses it,
never at file scope. A composition emits two copies of a body and namespaces
only the function declarations, so a file-scope const would be declared
twice by metamer x metamer and fail to link. */
/* One piecewise-continuous Gaussian lobe: the width differs either side of the
centre, which is how the fit follows an asymmetry a plain Gaussian cannot.
g and d are reciprocal widths. */
float meta_lobe(float lambda, float a, float b, float g, float d){
float x = (lambda - b) * (lambda < b ? g : d);
return a * exp(-0.5f * x * x);
}
/* x-bar, y-bar, z-bar of the CIE 1931 two-degree observer: three lobes, then
two and two. The negative amplitude is the real dip in x-bar near 500 nm. */
float3 meta_cmf(float lambda){
return float3(
meta_lobe(lambda, 0.362f, 442.0f, 0.0624f, 0.0374f)
+ meta_lobe(lambda, 1.056f, 599.8f, 0.0264f, 0.0323f)
+ meta_lobe(lambda, -0.065f, 501.1f, 0.0490f, 0.0382f),
meta_lobe(lambda, 0.821f, 568.8f, 0.0213f, 0.0247f)
+ meta_lobe(lambda, 0.286f, 530.9f, 0.0613f, 0.0322f),
meta_lobe(lambda, 1.217f, 437.0f, 0.0845f, 0.0278f)
+ meta_lobe(lambda, 0.681f, 459.0f, 0.0385f, 0.0725f));
}
/* CIE XYZ to linear sRGB (IEC 61966-2-1 primaries at D65), then that
standard's transfer function — a linear toe under a shifted power curve. */
float3 meta_rgb(float3 c){
return float3( 3.2404542f * c.x - 1.5371385f * c.y - 0.4985314f * c.z,
-0.9692660f * c.x + 1.8760108f * c.y + 0.0415560f * c.z,
0.0556434f * c.x - 0.2040259f * c.y + 1.0572252f * c.z);
}
float3 meta_encode(float3 c){
c = forma_clamp(c, 0.0f, 1.0f);
return lerp(12.92f * c, 1.055f * pow(c, float3(1.0f / 2.4f)) - 0.055f, forma_step(0.0031308f, c));
}
/* An exponential tilt across the band, peak normalised to 1.f Strictly
positive everywhere: the amplitude bound is set by wherever this is
smallest, and a base that touched zero would admit no black at all. */
float meta_base(float u, float tilt){
return exp(tilt * (u - 0.5f) - fabs(tilt) * 0.5f);
}
/* The discrete inner product of two vectors written in coefficient space over
(x-bar, y-bar, z-bar): everything the orthogonalisation needs to know about
the observer is the 3 x 3 Gram matrix of those three rows. */
float meta_ip(float3 a, float3 b, float3x3 G){ return dot(a, G * b); }
/* The metameric black at one wavelength. Five scalars are its whole
description — two coefficients on the test functions and a triple on the
matching functions — which is what lets a fragment shader evaluate it
anywhere without an array. The matching functions arrive already evaluated,
because every caller here has just needed them for something else and they
are seven exponentials each. */
float meta_black(float u, float3 c, float dir, float A1, float A2, float3 Q){
const float forma_PI = 3.141592653589793f;
return A1 * cos(dir * forma_PI * u) + A2 * cos((dir + 1.0f) * forma_PI * u) + dot(Q, c);
}
/* Chromaticity, not radiance: the base level of a spectrum is arbitrary here,
so each half of the band is scaled to the brightest sRGB triple carrying its
own colour. Measured across the whole declared box, nothing this plate can
mix is outside the gamut, so the forma_max() is a formality rather than a clamp. */
float3 meta_swatch(float3 xyz){
float3 lin = meta_rgb(xyz / forma_max(xyz.y, 1e-20));
return meta_encode(forma_max(lin, 0.0f) / forma_max(forma_max(lin.x, lin.y), forma_max(lin.z, 1e-9)));
}
/* 1 inside a mark and 0 outside it, with the transition spread over the
outer part of its width so the edge is antialiased rather than stepped. */
float meta_mark(float e0, float e1, float x){ return 1.0f - forma_smoothstep(e0, e1, x); }
float3 plate(float2 uv){
float3 forma_r = float3(0.0f, 0.0f, 0.0f);
for (int forma_once = 0; forma_once < 1; forma_once++){
/* The plate in three registers, as fractions of its height. */
const float PLOT_T = 0.055f;
const float PLOT_B = 0.648f;
const float RULE_T = 0.674f;
const float RULE_B = 0.708f;
const float BAND_T = 0.756f;
const float BAND_B = 0.960f;
const float TICK_A = 0.722f;
const float TICK_B = 0.744f;
const float TICK_C = 0.972f;
const float TICK_D = 0.994f;
const float LO = 380.0f;
const float SPAN = 400.0f;
/* The largest sRGB channel the matching functions reach when taken straight
through the matrix — measured over the band, at 605.7f nm. */
const float HUE_PEAK = 2.503f;
const float forma_PI = 3.141592653589793f;
float asp = u_res.y / u_res.x;
/* Four is the floor the mathematics has: a 3 x N matrix has no null space
at all until N exceeds 3.f The slider starts at 12 for a separate reason
(see the header), so this max is never the binding one — it is here so
that neither path can be handed a count that makes the object vanish. */
int n = forma_max(4, int(p_samples + 0.5f));
float fn = float(n);
float dir = floor(p_nulldir + 0.5f);
/* ---- pass one: the observer's Gram matrix, the test functions against it,
and the base spectrum's own tristimulus. Constant loop bound at the
slider's own ceiling, break on the live value, mandelbrot fashion. */
float3x3 G = float3x3(0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f);
float3 d1 = float3(0.0f);
float3 d2 = float3(0.0f);
float3 Sxyz = float3(0.0f);
float3 ff = float3(0.0f);
for (int k = 0; k < 88; k++){
if (k >= n) break;
float u = (float(k) + 0.5f) / fn;
float3 c = meta_cmf(LO + u * SPAN);
float s = meta_base(u, p_tilt);
float g1 = cos(dir * forma_PI * u);
float g2 = cos((dir + 1.0f) * forma_PI * u);
{ float3 forma_col0 = c.x * c; G[0][0] += forma_col0.x; G[1][0] += forma_col0.y; G[2][0] += forma_col0.z; } { float3 forma_col1 = c.y * c; G[0][1] += forma_col1.x; G[1][1] += forma_col1.y; G[2][1] += forma_col1.z; } { float3 forma_col2 = c.z * c; G[0][2] += forma_col2.x; G[1][2] += forma_col2.y; G[2][2] += forma_col2.z; }
d1 += g1 * c; d2 += g2 * c;
ff += float3(g1 * g1, g1 * g2, g2 * g2);
Sxyz += s * c;
}
/* Gram–Schmidt on the three rows, in coefficient space rather than over the
N samples: a vector is three numbers and the inner product is G. Cohen and
Kappauf write the same projector as A(A'A)inv A'; this is that operator
built without an inversion. */
float3 e1 = float3(1.0f, 0.0f, 0.0f);
e1 = e1 / sqrt(meta_ip(e1, e1, G));
float3 v2 = float3(0.0f, 1.0f, 0.0f);
v2 = v2 - e1 * meta_ip(v2, e1, G);
float3 e2 = v2 / sqrt(meta_ip(v2, v2, G));
float3 v3 = float3(0.0f, 0.0f, 1.0f);
v3 = v3 - e1 * meta_ip(v3, e1, G);
v3 = v3 - e2 * meta_ip(v3, e2, G);
float3 e3 = v3 / sqrt(meta_ip(v3, v3, G));
/* Projecting the two test functions off that basis leaves two vectors
orthogonal to all three rows, which is what being in the null space means.
Their mutual inner products come from numbers already in hand. */
float3 a1 = float3(dot(d1, e1), dot(d1, e2), dot(d1, e3));
float3 a2 = float3(dot(d2, e1), dot(d2, e2), dot(d2, e3));
float3 q1 = a1.x * e1 + a1.y * e2 + a1.z * e3;
float3 q2 = a2.x * e1 + a2.y * e2 + a2.z * e3;
float b11 = ff.x - dot(a1, a1);
float b12 = ff.y - dot(a1, a2);
float b22 = ff.z - dot(a2, a2);
float n1 = sqrt(forma_max(b11, 1e-12));
/* Orthonormalising the second against the first makes the drift a rotation
of a circle rather than a wander round an oblique pair. The floor guards
the case where the sample count cannot resolve the two test functions
apart, which is why the sample floor is 12; measured, the smallest this
ratio reaches inside the declared box is 0.81. */
float n2 = forma_max(sqrt(forma_max(b22 - b12 * b12 / (n1 * n1), 0.0f)), 1e-4 * n1);
float mu = b12 / (n1 * n1);
/* The drift: a rotation inside the null space. cycle() carries this plate's
own PHASE, so no two cards sit at the same angle. */
float th = 6.28318530718f * cycle(u_t, 40.0f);
float A1 = cos(th) / n1 - sin(th) * mu / n2;
float A2 = sin(th) / n2;
float3 Q = -(A1 * q1 + A2 * q2);
/* ---- pass two: the bound, and the black's own tristimulus. That second
sum is the claim, computed rather than assumed — if the construction is
right it is zero, and the two halves of the band are built by adding and
subtracting it. */
float limit = 1e20;
float3 D = float3(0.0f);
for (int k = 0; k < 88; k++){
if (k >= n) break;
float u = (float(k) + 0.5f) / fn;
float3 c = meta_cmf(LO + u * SPAN);
float b = meta_black(u, c, dir, A1, A2, Q);
/* A spectrum cannot carry negative light, so the multiple that may be
added is capped by whichever sample runs out of base spectrum first. */
if (fabs(b) > 1e-9) limit = forma_min(limit, meta_base(u, p_tilt) / fabs(b));
D += b * c;
}
float scale = p_amp * limit;
float PMAX = 1.0f + p_amp; // the base peak is 1, so this is the ceiling
/* ---- the exhibit ---------------------------------------------------- */
if (uv.y > BAND_T && uv.y < BAND_B){
/* Each half from its own spectrum. Factoring the constant out of the sum
is arithmetic, not a shortcut: sum of (s +- scale.z).cbar really is
sum of s.cbar +- scale times sum of b.cbar, and that second sum is the
quantity just measured rather than one assumed. The join is at 0.5f and
nothing marks it. */
{ forma_r = meta_swatch(uv.x < 0.5f ? Sxyz + scale * D : Sxyz - scale * D); break; }
}
if (uv.y > RULE_T && uv.y < RULE_B){
/* The wavelength rule, so the horizontal reads as wavelength: the matching
functions straight through the sRGB matrix — the fit used only for what
it was fitted to, the absolute values — with negatives clipped, because
a single wavelength is outside the gamut nearly everywhere. */
float3 l = meta_rgb(meta_cmf(LO + uv.x * SPAN));
{ forma_r = meta_encode(forma_max(l, 0.0f) / HUE_PEAK); break; }
}
if ((uv.y > TICK_A && uv.y < TICK_B) || (uv.y > TICK_C && uv.y < TICK_D)){
/* Two short ticks, above the band and below it, aimed at the join. */
{ forma_r = float3(0.016f, 0.024f, 0.039f)
+ float3(0.55f, 0.58f, 0.66f) * meta_mark(0.0018f, 0.0034f, fabs(uv.x - 0.5f)); break; }
}
if (uv.y > PLOT_T - 0.02f && uv.y < RULE_T){
float3 col = float3(4.0f, 6.0f, 10.0f) / 255.0f;
/* The zero axis, so that a curve touching it is legible as touching. */
col = lerp(col, float3(0.28f, 0.30f, 0.36f), meta_mark(0.0030f, 0.0052f, fabs(uv.y - PLOT_B) * asp));
/* The two spectra between the outer samples, straight-line interpolated
between them. The object is the dots; the line between them is scenery,
because only the samples are constrained. */
float x = uv.x * fn - 0.5f;
if (x >= 0.0f && x <= fn - 1.0f){
float k0f = forma_clamp(floor(x), 0.0f, fn - 2.0f);
float fr = forma_clamp(x - k0f, 0.0f, 1.0f);
float u0 = (k0f + 0.5f) / fn;
float u1 = (k0f + 1.5f) / fn;
float s0 = meta_base(u0, p_tilt);
float s1 = meta_base(u1, p_tilt);
float b0 = meta_black(u0, meta_cmf(LO + u0 * SPAN), dir, A1, A2, Q);
float b1 = meta_black(u1, meta_cmf(LO + u1 * SPAN), dir, A1, A2, Q);
/* The three curves at the two ends of this segment, packed so that one
line of arithmetic serves all of them: x is the base plus the black,
y the base minus it, z the base itself. */
float3 y0 = float3(s0 + scale * b0, s0 - scale * b0, s0);
float3 y1 = float3(s1 + scale * b1, s1 - scale * b1, s1);
float3 c0 = PLOT_B - (y0 / PMAX) * (PLOT_B - PLOT_T);
float3 c1 = PLOT_B - (y1 / PMAX) * (PLOT_B - PLOT_T);
/* A line drawn by vertical distance alone thins wherever it climbs, so
the vertical offset is turned into a perpendicular one through the
segment's own screen slope. */
float3 slope = (c1 - c0) * fn * asp;
float3 per = rsqrt(1.0f + slope * slope);
float3 vc = lerp(c0, c1, fr);
/* Sample dots, on the nearest sample's own centre — the object being
drawn is those points, and the line between them is scenery. */
float kd = forma_clamp(floor(uv.x * fn), 0.0f, fn - 1.0f);
float ud = (kd + 0.5f) / fn;
float sd = meta_base(ud, p_tilt);
float bd = meta_black(ud, meta_cmf(LO + ud * SPAN), dir, A1, A2, Q);
float2 dv = PLOT_B - (float2(sd + scale * bd, sd - scale * bd) / PMAX) * (PLOT_B - PLOT_T);
float du = uv.x - ud;
float LW = 0.0075f; // line half width, in frame widths
float DR = forma_min(0.0135f, 0.45f / fn); // dot radius, same units
col = lerp(col, float3(0.34f, 0.36f, 0.42f),
meta_mark(LW * 0.55f, LW, fabs(uv.y - vc.z) * asp * per.z) * 0.55f);
col = lerp(col, float3(0.80f, 0.85f, 0.94f),
forma_max(meta_mark(LW * 0.55f, LW, fabs(uv.y - vc.y) * asp * per.y),
meta_mark(DR * 0.6f, DR, length(float2(du, (uv.y - dv.y) * asp)))));
/* The order accent, in the ramp's bright lobe: the ramp is black around
0.5f and these are sparse strokes on a dark ground. */
col = lerp(col, ramp(0.92f),
forma_max(meta_mark(LW * 0.55f, LW, fabs(uv.y - vc.x) * asp * per.x),
meta_mark(DR * 0.6f, DR, length(float2(du, (uv.y - dv.x) * asp)))));
}
{ forma_r = col; break; }
}
{ forma_r = float3(4.0f, 6.0f, 10.0f) / 255.0f; break; }
}
return forma_r;
}
void process(int2 pos){
// FORMA's uv runs y-down like its canvas; Nuke's rows run up
float2 uv = float2((float(pos.x) + 0.5f) / u_res.x, 1.0f - (float(pos.y) + 0.5f) / u_res.y);
float3 c = plate(uv);
dst() = float4(c.x, c.y, c.z, 1.0f);
}
};