"""cvdcheck.py: lightness reversals, CVD simulation and contrast for any legend. Single file, numpy only. Give it an ordered list of sRGB colours (low value to high value) and it reports: * CIE L* profile (sRGB -> linear -> XYZ (D65) -> CIELAB) * number of L* reversals: turning points whose swing exceeds a threshold (hysteresis count, so slow drifts in 256-entry tables are not lost) * CIEDE2000 between adjacent classes, under normal vision and under Machado, Oliveira & Fernandes (2009) severity-1.0 protan, deutan and tritan simulation (matrices applied to linear sRGB, then clipped) * WCAG 2 contrast ratio for adjacent pairs and for all pairs Usage: python cvdcheck.py --selftest python cvdcheck.py "#ff0000" "#00ff00" "#0000ff" Lea Keller, agentik.blog Lab, 2026-10-04. """ import sys import numpy as np # ---------------------------------------------------------------- colour M_RGB2XYZ = np.array([[0.4124564, 0.3575761, 0.1804375], [0.2126729, 0.7151522, 0.0721750], [0.0193339, 0.1191920, 0.9503041]]) WHITE = M_RGB2XYZ @ np.ones(3) # D65 white of the sRGB primaries # Machado et al. 2009, severity 1.0 (supplementary tables, as in colorspacious) MACHADO = { "protan": np.array([[0.152286, 1.052583, -0.204868], [0.114503, 0.786281, 0.099216], [-0.003882, -0.048116, 1.051998]]), "deutan": np.array([[0.367322, 0.860646, -0.227968], [0.280085, 0.672501, 0.047413], [-0.011820, 0.042940, 0.968881]]), "tritan": np.array([[1.255528, -0.076749, -0.178779], [-0.078411, 0.930809, 0.147602], [0.004733, 0.691367, 0.303900]]), } def parse_color(c): """'#rrggbb', (r,g,b) in 0..1, or (r,g,b) in 0..255 -> float array 0..1.""" if isinstance(c, str): c = c.strip().strip('"').lstrip('#') return np.array([int(c[i:i + 2], 16) for i in (0, 2, 4)]) / 255.0 c = np.asarray(c, float) return c / 255.0 if c.max() > 1.0 else c def srgb_to_linear(rgb): rgb = np.asarray(rgb, float) return np.where(rgb <= 0.04045, rgb / 12.92, ((rgb + 0.055) / 1.055) ** 2.4) def linear_to_srgb(lin): lin = np.clip(lin, 0, 1) return np.where(lin <= 0.0031308, 12.92 * lin, 1.055 * lin ** (1 / 2.4) - 0.055) def linear_to_lab(lin): xyz = lin @ M_RGB2XYZ.T / WHITE d = 6 / 29 f = np.where(xyz > d ** 3, np.cbrt(xyz), xyz / (3 * d * d) + 4 / 29) L = 116 * f[..., 1] - 16 a = 500 * (f[..., 0] - f[..., 1]) b = 200 * (f[..., 1] - f[..., 2]) return np.stack([L, a, b], -1) def srgb_to_lab(rgb): return linear_to_lab(srgb_to_linear(rgb)) def simulate(rgb, kind): """Machado 2009 severity 1.0 in linear RGB; returns sRGB 0..1.""" if kind == "normal": return np.asarray(rgb, float) lin = srgb_to_linear(rgb) @ MACHADO[kind].T return linear_to_srgb(np.clip(lin, 0, 1)) # ---------------------------------------------------------------- CIEDE2000 def ciede2000(lab1, lab2): """Sharma, Wu & Dalal (2005) implementation, kL = kC = kH = 1.""" lab1, lab2 = np.atleast_2d(lab1).astype(float), np.atleast_2d(lab2).astype(float) L1, a1, b1 = lab1.T L2, a2, b2 = lab2.T C1, C2 = np.hypot(a1, b1), np.hypot(a2, b2) Cb = (C1 + C2) / 2 G = 0.5 * (1 - np.sqrt(Cb ** 7 / (Cb ** 7 + 25.0 ** 7))) a1p, a2p = (1 + G) * a1, (1 + G) * a2 C1p, C2p = np.hypot(a1p, b1), np.hypot(a2p, b2) h1p = np.degrees(np.arctan2(b1, a1p)) % 360 h2p = np.degrees(np.arctan2(b2, a2p)) % 360 h1p = np.where((b1 == 0) & (a1p == 0), 0, h1p) h2p = np.where((b2 == 0) & (a2p == 0), 0, h2p) dLp = L2 - L1 dCp = C2p - C1p dh = h2p - h1p dh = np.where(dh > 180, dh - 360, dh) dh = np.where(dh < -180, dh + 360, dh) dh = np.where(C1p * C2p == 0, 0, dh) dHp = 2 * np.sqrt(C1p * C2p) * np.sin(np.radians(dh / 2)) Lbp = (L1 + L2) / 2 Cbp = (C1p + C2p) / 2 hs = h1p + h2p hbp = np.where(np.abs(h1p - h2p) > 180, np.where(hs < 360, (hs + 360) / 2, (hs - 360) / 2), hs / 2) hbp = np.where(C1p * C2p == 0, hs, hbp) T = (1 - 0.17 * np.cos(np.radians(hbp - 30)) + 0.24 * np.cos(np.radians(2 * hbp)) + 0.32 * np.cos(np.radians(3 * hbp + 6)) - 0.20 * np.cos(np.radians(4 * hbp - 63))) dtheta = 30 * np.exp(-((hbp - 275) / 25) ** 2) Rc = 2 * np.sqrt(Cbp ** 7 / (Cbp ** 7 + 25.0 ** 7)) Sl = 1 + 0.015 * (Lbp - 50) ** 2 / np.sqrt(20 + (Lbp - 50) ** 2) Sc = 1 + 0.045 * Cbp Sh = 1 + 0.015 * Cbp * T Rt = -np.sin(np.radians(2 * dtheta)) * Rc return np.sqrt((dLp / Sl) ** 2 + (dCp / Sc) ** 2 + (dHp / Sh) ** 2 + Rt * (dCp / Sc) * (dHp / Sh)) # ---------------------------------------------------------------- WCAG 2 def rel_luminance(rgb): # WCAG 2.x relative luminance (2.2 text uses the 0.04045 knee, as here) return srgb_to_linear(rgb) @ np.array([0.2126, 0.7152, 0.0722]) def wcag_ratio(rgb1, rgb2): l1, l2 = rel_luminance(rgb1), rel_luminance(rgb2) hi, lo = np.maximum(l1, l2), np.minimum(l1, l2) return (hi + 0.05) / (lo + 0.05) # ---------------------------------------------------------------- reversals def count_reversals(L, thr=1.0): """Turning points in a sequence whose swing is at least thr L* units. Hysteresis: a reversal is counted when L moves thr or more against the current direction, measured from the most recent extreme. Plateaus and wiggles smaller than thr are ignored. Returns (count, indices of extremes). """ L = np.asarray(L, float) direction, ext_i, n, turns = 0, 0, 0, [] for i in range(1, len(L)): d = L[i] - L[ext_i] if direction == 0: if abs(d) >= thr: direction = np.sign(d) ext_i = i elif direction > 0: if L[i] > L[ext_i]: ext_i = i elif L[ext_i] - L[i] >= thr: n += 1; turns.append(ext_i); direction = -1; ext_i = i else: if L[i] < L[ext_i]: ext_i = i elif L[i] - L[ext_i] >= thr: n += 1; turns.append(ext_i); direction = 1; ext_i = i return n, turns # ---------------------------------------------------------------- report VISIONS = ("normal", "deutan", "protan", "tritan") def analyse(colors, de_floor=5.0): rgb = np.array([parse_color(c) for c in colors]) out = {"n": len(rgb)} lab = srgb_to_lab(rgb) out["L"] = lab[:, 0] for t in (0.5, 1.0, 2.0): out[f"rev_{t}"] = count_reversals(lab[:, 0], t)[0] for v in VISIONS: s = simulate(rgb, v) sl = srgb_to_lab(s) de = ciede2000(sl[:-1], sl[1:]) out[f"de_{v}"] = de out[f"rev1_{v}"] = count_reversals(sl[:, 0], 1.0)[0] out[f"min_de_{v}"] = de.min() out[f"n_de_under_{v}"] = int((de < de_floor).sum()) out["collapse_deutan"] = int(((out["de_normal"] > 10) & (out["de_deutan"] < 5)).sum()) out["collapse_protan"] = int(((out["de_normal"] > 10) & (out["de_protan"] < 5)).sum()) adj = wcag_ratio(rgb[:-1], rgb[1:]) out["wcag_adj"] = adj out["min_wcag_adj"] = adj.min() i, j = np.triu_indices(len(rgb), 1) allp = wcag_ratio(rgb[i], rgb[j]) out["min_wcag_all"] = allp.min() out["frac_all_pairs_wcag_under_1.5"] = float((allp < 1.5).mean()) return out # ---------------------------------------------------------------- self-test SHARMA_SOURCE = "data/gfiumara_testCIEDE2000.cpp" def selftest(sharma_path=SHARMA_SOURCE, colorspacious_path="data/colorspacious_cvd.py"): import re ok = True r = wcag_ratio(np.array([1., 1, 1]), np.array([0., 0, 0])) print(f"WCAG white/black = {r:.6f} (expect 21)", "PASS" if abs(r - 21) < 1e-4 else "FAIL") ok &= abs(r - 21) < 1e-4 lw = srgb_to_lab(np.array([1., 1, 1])) print(f"Lab(white) = {np.round(lw, 6)} (expect 100,0,0)", "PASS" if np.allclose(lw, [100, 0, 0], atol=1e-4) else "FAIL") ok &= np.allclose(lw, [100, 0, 0], atol=1e-4) # Machado: compare with an independent copy of the supplementary tables txt = open(colorspacious_path).read() names = {"protan": "protanomaly", "deutan": "deuteranomaly", "tritan": "tritanomaly"} for k, name in names.items(): dict_start = txt.index("MACHADO_ET_AL_MATRICES = {") block = txt[txt.index(f'"{name}"', dict_start):] m = re.search(r"100:\s*\[(.*?)\]\s*,\s*\]", block, re.S) vals = np.array([float(x) for x in re.findall(r"-?\d+\.\d+", m.group(1))]).reshape(3, 3) good = np.allclose(vals, MACHADO[k], atol=1e-6) print(f"Machado {k} 1.0 matrix vs colorspacious copy: max diff " f"{np.abs(vals - MACHADO[k]).max():.1e}", "PASS" if good else "FAIL") ok &= good red = srgb_to_linear(np.array([1., 0, 0])) @ MACHADO["deutan"].T good = np.allclose(red, [0.367322, 0.280085, -0.011820], atol=1e-6) print(f"deutan(pure red), linear = {np.round(red, 6)} (expect first column)", "PASS" if good else "FAIL") ok &= good # CIEDE2000: 34 pairs of Sharma, Wu & Dalal (2005), via gfiumara/CIEDE2000 src = open(sharma_path).read() l1 = re.findall(r"lab1 = \{([^}]*)\}", src) l2 = re.findall(r"lab2 = \{([^}]*)\}", src) ex = re.findall(r"expectedResult = ([\d.]+);", src) A = np.array([[float(x) for x in s.split(",")] for s in l1]) B = np.array([[float(x) for x in s.split(",")] for s in l2]) E = np.array([float(x) for x in ex]) D = ciede2000(A, B) D2 = ciede2000(B, A) err = np.abs(D - E).max() sym = np.abs(D - D2).max() good = (len(E) == 34) and err <= 1e-4 and sym < 1e-9 print(f"CIEDE2000: {len(E)} Sharma pairs, max |error| = {err:.2e}, " f"max asymmetry = {sym:.1e}", "PASS" if good else "FAIL") ok &= good print("ALL PASS" if ok else "SELFTEST FAILED") return ok if __name__ == "__main__": if len(sys.argv) > 1 and sys.argv[1] == "--selftest": sys.exit(0 if selftest() else 1) res = analyse(sys.argv[1:]) print(f"classes {res['n']}; L* reversals (0.5/1/2): " f"{res['rev_0.5']}/{res['rev_1.0']}/{res['rev_2.0']}") print("L*:", np.round(res["L"], 1)) for v in VISIONS: print(f"{v:7s} adjacent dE2000:", np.round(res[f"de_{v}"], 1)) print("adjacent WCAG ratios:", np.round(res["wcag_adj"], 2))