nx_colorsci_lib.nx source
↩ module page · 872 lines · 38414 B
1// nx_colorsci_lib.nx -- COLOUR SCIENCE AS ONE SOVEREIGN RULER (2026-08-24, phototwin PT1: cs_de2000, cs_icc_gamut).
2// license_tier: ORIGINAL No hw writes (Rule 26). LIB (no main).
3//
4// WHAT IT IS. CIE 1976 L*a*b* both ways (CIE 15:2004), CIEDE2000 (CIE 142-2001, Sharma-Wu-Dalal 2005 eqs 4-22,
5// KAT = their Table I as data in knowledge/colorsci_ciede2000_kat.tsv), an ICC v2/v4 matrix/TRC profile reader
6// (rXYZ gXYZ bXYZ wtpt rTRC gTRC bTRC; the PCS is D50 by the spec, so image and substrate meet in ONE space), a
7// segment-maxima gamut boundary descriptor (Morovic and Luo 2000) and chroma clipping toward the achromatic
8// axis at constant L and hue. All arithmetic is nx_fixq30_lib Q30 fixed point: no float, no typed-in constants.
9//
10// WHY THROUGH THE PROFILE. The estate had RGB/HSV/YCbCr and no CIELAB, so no colour difference could be stated in
11// the unit the whole field uses. Going RGB -> TRC -> matrix -> XYZ(D50) -> Lab through a parsed ICC profile is
12// the ICC workflow itself, which means a measured substrate profile from a swatch bench drops in with no code.
13//
14// EVERY THRESHOLD IS CONF (knowledge/colorsci.conf, cs_conf_int) and an absent key is a REFUSAL, never a default.
15import "nx_syscalls.nx"
16import "nx_fixq30_lib.nx"
17
18// ---- CIE 1976 L*a*b* ----
19const CS_LAB_DELTA_NUM: i64 = 6
20const CS_LAB_DELTA_DEN: i64 = 29
21const CS_LAB_L_SCALE: i64 = 116
22const CS_LAB_L_OFFSET: i64 = 16
23const CS_LAB_A_SCALE: i64 = 500
24const CS_LAB_B_SCALE: i64 = 200
25const CS_LAB_F_LIN_NUM: i64 = 4
26const CS_LAB_F_LIN_DEN: i64 = 29
27const CS_THREE: i64 = 3
28const CS_TWO: i64 = 2
29// ---- CIEDE2000 ----
30const CS_DE_C_REF: i64 = 25
31const CS_DE_T_W1_MILLI: i64 = 170
32const CS_DE_T_W2_MILLI: i64 = 240
33const CS_DE_T_W3_MILLI: i64 = 320
34const CS_DE_T_W4_MILLI: i64 = 200
35const CS_DE_T_D1_DEG: i64 = 30
36const CS_DE_T_M2: i64 = 2
37const CS_DE_T_M3: i64 = 3
38const CS_DE_T_D3_DEG: i64 = 6
39const CS_DE_T_M4: i64 = 4
40const CS_DE_T_D4_DEG: i64 = 63
41const CS_DE_THETA_CENTER_DEG: i64 = 275
42const CS_DE_THETA_WIDTH_DEG: i64 = 25
43const CS_DE_THETA_SCALE_DEG: i64 = 30
44const CS_DE_SL_K_MILLI: i64 = 15
45const CS_DE_SL_L0: i64 = 50
46const CS_DE_SL_C: i64 = 20
47const CS_DE_SC_K_MILLI: i64 = 45
48const CS_DE_SH_K_MILLI: i64 = 15
49const CS_MILLI: i64 = 1000
50const CS_DEG_FULL: i64 = 360
51const CS_DEG_HALF: i64 = 180
52const CS_HUE_KNIFE_EPS_Q30: i64 = 64 // 6e-8 degrees: the rounding bound of two Q30 rad->deg conversions (see cs_de2000_ex)
53// ---- context slots ----
54const CS_X_FQ: i64 = 0
55const CS_X_DELTA3: i64 = 1
56const CS_X_INV3D2: i64 = 2
57const CS_X_FLIN_OFF: i64 = 3
58const CS_X_DELTA: i64 = 4
59const CS_X_SLOTS: i64 = 8
60const CS_SLOT: i64 = 8
61// ---- ICC (ICC.1:2010) ----
62const CS_ICC_HDR_BYTES: i64 = 128
63const CS_ICC_MAGIC_OFF: i64 = 36
64const CS_ICC_TAG_ENTRY_BYTES: i64 = 12
65const CS_ICC_TAG_OFF_FIELD: i64 = 4
66const CS_ICC_TAG_SIZE_FIELD: i64 = 8
67const CS_ICC_TYPE_HDR_BYTES: i64 = 8
68const CS_ICC_XYZ_BYTES: i64 = 4
69const CS_ICC_S15F16_SCALE: i64 = 16384 // 2^(30-16): s15Fixed16 -> Q30
70const CS_ICC_U8F8_SCALE: i64 = 4194304 // 2^(30-8): u8Fixed8 -> Q30
71const CS_ICC_U16_MAX: i64 = 65535
72const CS_ICC_U32_SIGN: i64 = 2147483648
73const CS_ICC_U32_WRAP: i64 = 4294967296
74const CS_ICC_TRC_IDENTITY: i64 = 0
75const CS_ICC_TRC_GAMMA: i64 = 1
76const CS_ICC_TRC_TABLE: i64 = 2
77const CS_ICC_TAGS_NEEDED: i64 = 7
78const CS_BYTE_BITS_3: i64 = 24
79const CS_BYTE_BITS_2: i64 = 16
80const CS_BYTE_BITS_1: i64 = 8
81const CS_BYTE_MAX: i64 = 255
82const CS_LUT_ENTRIES: i64 = 256
83const CS_FOURCC: i64 = 4
84const CS_P_MAT: i64 = 0
85const CS_P_WTPT: i64 = 9
86const CS_P_LUT_R: i64 = 12
87const CS_P_LUT_G: i64 = 13
88const CS_P_LUT_B: i64 = 14
89const CS_P_INV: i64 = 15
90const CS_P_OK: i64 = 24
91const CS_P_BYTES: i64 = 25
92const CS_P_TAGS: i64 = 26
93const CS_P_TRC_R: i64 = 27
94const CS_P_TRC_G: i64 = 28
95const CS_P_TRC_B: i64 = 29
96const CS_P_HAVE: i64 = 30
97const CS_P_B64: i64 = 31 // 1 when the profile arrived as base64 text
98const CS_P_SLOTS: i64 = 32
99// base64 alphabet (RFC 4648 section 4)
100const CS_B64_UPPER_A: i64 = 65
101const CS_B64_UPPER_Z: i64 = 90
102const CS_B64_LOWER_A: i64 = 97
103const CS_B64_LOWER_Z: i64 = 122
104const CS_B64_LOWER_BASE: i64 = 26
105const CS_B64_DIGIT_BASE: i64 = 52
106const CS_B64_PLUS: i64 = 43
107const CS_B64_PLUS_VAL: i64 = 62
108const CS_B64_SLASH: i64 = 47
109const CS_B64_SLASH_VAL: i64 = 63
110const CS_B64_PAD: i64 = 61
111const CS_B64_BITS: i64 = 6
112const CS_MAT_N: i64 = 9
113const CS_MAT_COLS: i64 = 3
114// ---- gamut boundary descriptor ----
115const CS_G_H: i64 = 0
116const CS_G_E: i64 = 1
117const CS_G_COUNT: i64 = 2
118const CS_G_BASE: i64 = 3
119const CS_G_CENTER_L: i64 = 50
120const CS_PERMIL: i64 = 1000
121const CS_G_UNSAMPLED: i64 = 0 - 1
122const CS_CUBE_FACES: i64 = 6
123const CS_CLIP_LOW_RADIUS: i64 = 0
124// ---- conf reader ----
125const CS_CONF_MISS: i64 = 0 - 999999999
126const CS_ASCII_NL: i64 = 10
127const CS_ASCII_EQ: i64 = 61
128const CS_ASCII_ZERO: i64 = 48
129const CS_ASCII_NINE: i64 = 57
130const CS_DECIMAL: i64 = 10
131
132// ---- CHROMATIC ADAPTATION (CIE 15:2004 Annex E; Bradford is the transform the ICC PCS is defined against,
133// CAT02 is its CIECAM02 successor). WHY HERE AND NOT IN A SECOND LIB: colour science is ONE ruler in this
134// estate, and an adaptation living beside the profile reader is the only way a D65 capture and the D50 PCS
135// can meet without a second, drifting copy of the same arithmetic.
136// WHAT IT FIXES: the profile path applies matrix plus TRC only, so a colour measured under a non-PCS
137// illuminant carries a white-point error that no gamut work downstream can undo. The error is invisible
138// because every number stays in range and merely reads slightly wrong.
139// Coefficients are the published matrices carried at 1/10000, the same idiom as the CIEDE2000 weights above.
140const CS_CAT_SCALE: i64 = 10000
141const CS_CAT_BRADFORD: i64 = 0
142const CS_CAT_CAT02: i64 = 1
143const CS_CAT_M: i64 = 9
144const CS_CAT_BFD_0: i64 = 8951
145const CS_CAT_BFD_1: i64 = 2664
146const CS_CAT_BFD_2: i64 = 0 - 1614
147const CS_CAT_BFD_3: i64 = 0 - 7502
148const CS_CAT_BFD_4: i64 = 17135
149const CS_CAT_BFD_5: i64 = 367
150const CS_CAT_BFD_6: i64 = 389
151const CS_CAT_BFD_7: i64 = 0 - 685
152const CS_CAT_BFD_8: i64 = 10296
153const CS_CAT_02_0: i64 = 7328
154const CS_CAT_02_1: i64 = 4296
155const CS_CAT_02_2: i64 = 0 - 1624
156const CS_CAT_02_3: i64 = 0 - 7036
157const CS_CAT_02_4: i64 = 16975
158const CS_CAT_02_5: i64 = 61
159const CS_CAT_02_6: i64 = 30
160const CS_CAT_02_7: i64 = 136
161const CS_CAT_02_8: i64 = 9834
162
163// ---- byte helpers ----
164func cs_be32(b: *u8, off: i64) -> i64 {
165 return ((b[off] as i64) << CS_BYTE_BITS_3) + ((b[off + 1] as i64) << CS_BYTE_BITS_2) + ((b[off + 2] as i64) << CS_BYTE_BITS_1) + (b[off + 3] as i64)
166}
167func cs_be16(b: *u8, off: i64) -> i64 { return ((b[off] as i64) << CS_BYTE_BITS_1) + (b[off + 1] as i64) }
168func cs_s15f16_q30(b: *u8, off: i64) -> i64 {
169 var v: i64 = cs_be32(b, off)
170 if v >= CS_ICC_U32_SIGN { v = v - CS_ICC_U32_WRAP }
171 return v * CS_ICC_S15F16_SCALE
172}
173func cs_tag_is(b: *u8, off: i64, s: *u8) -> i64 {
174 var i: i64 = 0
175 while i < CS_FOURCC { if b[off + i] != s[i] { return 0 } i = i + 1 }
176 return 1
177}
178
179// ---- line-anchored conf integer read: CS_CONF_MISS when the key is absent (a comment cannot satisfy it) ----
180func cs_conf_int(path: *u8, key: *u8) -> i64 {
181 let lp: *i64 = sys_mmap(CS_SLOT * CS_TWO) as *i64
182 let b: *u8 = sys_read_file(path, lp)
183 if (b as i64) == 0 { return CS_CONF_MISS }
184 let n: i64 = lp[0]
185 var i: i64 = 0
186 var bol: i64 = 1
187 var result: i64 = CS_CONF_MISS
188 while i < n {
189 if result == CS_CONF_MISS { if bol == 1 {
190 var k: i64 = 0
191 var same: i64 = 1
192 while key[k] != (0 as u8) {
193 if i + k >= n { same = 0 }
194 if same == 1 { if b[i + k] != key[k] { same = 0 } }
195 k = k + 1
196 }
197 if same == 1 { if i + k < n { if (b[i + k] as i64) == CS_ASCII_EQ {
198 var o: i64 = i + k + 1
199 var v: i64 = 0
200 var digits: i64 = 0
201 var go: i64 = 1
202 while go == 1 {
203 if o >= n { go = 0 } else {
204 let c: i64 = b[o] as i64
205 if c < CS_ASCII_ZERO { go = 0 }
206 if c > CS_ASCII_NINE { go = 0 }
207 if go == 1 { v = v * CS_DECIMAL + (c - CS_ASCII_ZERO); digits = digits + 1; o = o + 1 }
208 }
209 }
210 if digits > 0 { result = v }
211 } } }
212 } }
213 if (b[i] as i64) == CS_ASCII_NL { bol = 1 } else { bol = 0 }
214 i = i + 1
215 }
216 sys_free_file(b, n)
217 return result
218}
219
220// ---- context ----
221func cs_ctx() -> *i64 {
222 let c: *i64 = sys_mmap(CS_X_SLOTS * CS_SLOT) as *i64
223 let fq: *i64 = fq_ctx()
224 c[CS_X_FQ] = fq as i64
225 let delta: i64 = fq_div(fq_from_int(CS_LAB_DELTA_NUM), fq_from_int(CS_LAB_DELTA_DEN))
226 c[CS_X_DELTA] = delta
227 let d2: i64 = fq_mul(delta, delta)
228 c[CS_X_DELTA3] = fq_mul(d2, delta)
229 c[CS_X_INV3D2] = fq_div(FQ_ONE, CS_THREE * d2)
230 c[CS_X_FLIN_OFF] = fq_div(fq_from_int(CS_LAB_F_LIN_NUM), fq_from_int(CS_LAB_F_LIN_DEN))
231 return c
232}
233func cs_fq(cs: *i64) -> *i64 { return cs[CS_X_FQ] as *i64 }
234
235// ---- CIELAB companding ----
236func cs_lab_f(cs: *i64, t: i64) -> i64 {
237 if t > cs[CS_X_DELTA3] { return fq_cbrt(cs_fq(cs), t) }
238 return fq_mul(t, cs[CS_X_INV3D2]) + cs[CS_X_FLIN_OFF]
239}
240func cs_lab_finv(cs: *i64, f: i64) -> i64 {
241 if f > cs[CS_X_DELTA] { return fq_mul(fq_mul(f, f), f) }
242 return fq_div(f - cs[CS_X_FLIN_OFF], cs[CS_X_INV3D2])
243}
244// XYZ (Q30, absolute) + white (Q30) -> L a b (Q30) in out3
245func cs_lab_from_xyz(cs: *i64, X: i64, Y: i64, Z: i64, Xn: i64, Yn: i64, Zn: i64, out3: *i64) -> i64 {
246 let fx: i64 = cs_lab_f(cs, fq_div(X, Xn))
247 let fy: i64 = cs_lab_f(cs, fq_div(Y, Yn))
248 let fz: i64 = cs_lab_f(cs, fq_div(Z, Zn))
249 out3[0] = CS_LAB_L_SCALE * fy - fq_from_int(CS_LAB_L_OFFSET)
250 out3[1] = CS_LAB_A_SCALE * (fx - fy)
251 out3[2] = CS_LAB_B_SCALE * (fy - fz)
252 return 0
253}
254func cs_xyz_from_lab(cs: *i64, L: i64, a: i64, b: i64, Xn: i64, Yn: i64, Zn: i64, out3: *i64) -> i64 {
255 let fy: i64 = fq_div(L + fq_from_int(CS_LAB_L_OFFSET), fq_from_int(CS_LAB_L_SCALE))
256 let fx: i64 = fy + fq_div(a, fq_from_int(CS_LAB_A_SCALE))
257 let fz: i64 = fy - fq_div(b, fq_from_int(CS_LAB_B_SCALE))
258 out3[0] = fq_mul(cs_lab_finv(cs, fx), Xn)
259 out3[1] = fq_mul(cs_lab_finv(cs, fy), Yn)
260 out3[2] = fq_mul(cs_lab_finv(cs, fz), Zn)
261 return 0
262}
263
264// ---- CIEDE2000 ----
265// sqrt(C^7 / (C^7 + 25^7)) with u = C/25, computed as (u^3 * sqrt(u)) / sqrt(u^7 + 1). MEASURED 2026-08-24 on the KAT:
266// forming u^7 first and taking the sqrt afterwards UNDERFLOWS Q30 for small chroma (C = 1.12 gives u^7 = 3.6e-10,
267// below one Q30 unit, so G read exactly 0.5 and the hue of Table I pairs 7/8 came out 1.7e-4 degrees high). The
268// numerator's factors are each representable; u <= 8 keeps u^7 under 2^51 so the denominator never overflows.
269func cs_de_c7_sqrt_ratio(fq: *i64, C: i64) -> i64 {
270 let u: i64 = fq_div(C, fq_from_int(CS_DE_C_REF))
271 if u == 0 { return 0 }
272 let u2: i64 = fq_mul(u, u)
273 let u3: i64 = fq_mul(u2, u)
274 let u4: i64 = fq_mul(u2, u2)
275 let u7: i64 = fq_mul(u4, u3)
276 let num: i64 = fq_mul(u3, fq_sqrt(u))
277 let den: i64 = fq_sqrt(u7 + FQ_ONE)
278 return fq_div(num, den)
279}
280// hue angle in DEGREES (Q30) in [0, 360); 0 when a = b = 0
281func cs_hue_deg(fq: *i64, a: i64, b: i64) -> i64 {
282 if a == 0 { if b == 0 { return 0 } }
283 var h: i64 = fq_rad2deg(fq, fq_atan2(fq, b, a))
284 if h < 0 { h = h + fq_from_int(CS_DEG_FULL) }
285 if h >= fq_from_int(CS_DEG_FULL) { h = h - fq_from_int(CS_DEG_FULL) }
286 return h
287}
288func cs_cos_deg(fq: *i64, deg: i64) -> i64 { return fq_cos(fq, fq_deg2rad(fq, deg)) }
289func cs_sin_deg(fq: *i64, deg: i64) -> i64 { return fq_sin(fq, fq_deg2rad(fq, deg)) }
290func cs_milli(m: i64) -> i64 { return fq_div(fq_from_int(m), fq_from_int(CS_MILLI)) }
291// THE RULER. All inputs Q30 in real CIELAB units. out4: [0] dE00, [1] h_bar_prime (degrees), [2] G, [3] T.
292func cs_de2000_ex(cs: *i64, L1: i64, a1: i64, b1: i64, L2: i64, a2: i64, b2: i64, out4: *i64) -> i64 {
293 let fq: *i64 = cs_fq(cs)
294 let d180: i64 = fq_from_int(CS_DEG_HALF)
295 let d360: i64 = fq_from_int(CS_DEG_FULL)
296 let C1: i64 = fq_hypot(a1, b1)
297 let C2: i64 = fq_hypot(a2, b2)
298 let Cbar: i64 = (C1 + C2) / CS_TWO
299 let G: i64 = (FQ_ONE - cs_de_c7_sqrt_ratio(fq, Cbar)) / CS_TWO
300 let a1p: i64 = fq_mul(FQ_ONE + G, a1)
301 let a2p: i64 = fq_mul(FQ_ONE + G, a2)
302 let C1p: i64 = fq_hypot(a1p, b1)
303 let C2p: i64 = fq_hypot(a2p, b2)
304 let h1p: i64 = cs_hue_deg(fq, a1p, b1)
305 let h2p: i64 = cs_hue_deg(fq, a2p, b2)
306 let dLp: i64 = L2 - L1
307 let dCp: i64 = C2p - C1p
308 var dhp: i64 = 0
309 var both: i64 = 0
310 if C1p != 0 { if C2p != 0 { both = 1 } }
311 // THE 180-DEGREE KNIFE EDGE (Table I pairs 9-16 exist to trap it): Sharma's branch is on |h1'-h2'| <= 180 EXACTLY,
312 // and pairs 10 and 14 sit on it. In Q30 the two hues round independently by a unit or two, so the raw compare
313 // took the wrong branch and returned 7.2195 for a published 7.1792. The named epsilon below is the bound of that
314 // rounding (two rad->deg conversions, a few units each), 6e-8 degrees -- six orders below the table's resolution.
315 let knife: i64 = d180 + CS_HUE_KNIFE_EPS_Q30
316 if both == 1 {
317 dhp = h2p - h1p
318 if dhp > knife { dhp = dhp - d360 }
319 if dhp < (0 - knife) { dhp = dhp + d360 }
320 }
321 let dHp: i64 = CS_TWO * fq_mul(fq_sqrt(fq_mul(C1p, C2p)), cs_sin_deg(fq, dhp / CS_TWO))
322 let Lbp: i64 = (L1 + L2) / CS_TWO
323 let Cbp: i64 = (C1p + C2p) / CS_TWO
324 var hbp: i64 = h1p + h2p
325 if both == 1 {
326 if fq_abs(h1p - h2p) <= knife { hbp = (h1p + h2p) / CS_TWO }
327 else { if h1p + h2p < d360 { hbp = (h1p + h2p + d360) / CS_TWO } else { hbp = (h1p + h2p - d360) / CS_TWO } }
328 }
329 let t1: i64 = fq_mul(cs_milli(CS_DE_T_W1_MILLI), cs_cos_deg(fq, hbp - fq_from_int(CS_DE_T_D1_DEG)))
330 let t2: i64 = fq_mul(cs_milli(CS_DE_T_W2_MILLI), cs_cos_deg(fq, CS_DE_T_M2 * hbp))
331 let t3: i64 = fq_mul(cs_milli(CS_DE_T_W3_MILLI), cs_cos_deg(fq, CS_DE_T_M3 * hbp + fq_from_int(CS_DE_T_D3_DEG)))
332 let t4: i64 = fq_mul(cs_milli(CS_DE_T_W4_MILLI), cs_cos_deg(fq, CS_DE_T_M4 * hbp - fq_from_int(CS_DE_T_D4_DEG)))
333 let T: i64 = FQ_ONE - t1 + t2 + t3 - t4
334 let yth: i64 = fq_div(hbp - fq_from_int(CS_DE_THETA_CENTER_DEG), fq_from_int(CS_DE_THETA_WIDTH_DEG))
335 let dtheta: i64 = CS_DE_THETA_SCALE_DEG * fq_exp(fq, 0 - fq_mul(yth, yth))
336 let RC: i64 = CS_TWO * cs_de_c7_sqrt_ratio(fq, Cbp)
337 let dL50: i64 = Lbp - fq_from_int(CS_DE_SL_L0)
338 let dL50sq: i64 = fq_mul(dL50, dL50)
339 let SL: i64 = FQ_ONE + fq_div(fq_mul(cs_milli(CS_DE_SL_K_MILLI), dL50sq), fq_sqrt(fq_from_int(CS_DE_SL_C) + dL50sq))
340 let SC: i64 = FQ_ONE + fq_mul(cs_milli(CS_DE_SC_K_MILLI), Cbp)
341 let SH: i64 = FQ_ONE + fq_mul(fq_mul(cs_milli(CS_DE_SH_K_MILLI), Cbp), T)
342 let RT: i64 = 0 - fq_mul(cs_sin_deg(fq, CS_TWO * dtheta), RC)
343 let tL: i64 = fq_div(dLp, SL)
344 let tC: i64 = fq_div(dCp, SC)
345 let tH: i64 = fq_div(dHp, SH)
346 let sum: i64 = fq_mul(tL, tL) + fq_mul(tC, tC) + fq_mul(tH, tH) + fq_mul(RT, fq_mul(tC, tH))
347 out4[0] = fq_sqrt(sum)
348 out4[1] = hbp
349 out4[2] = G
350 out4[3] = T
351 return 0
352}
353func cs_de2000(cs: *i64, L1: i64, a1: i64, b1: i64, L2: i64, a2: i64, b2: i64) -> i64 {
354 let o: *i64 = sys_mmap(CS_SLOT * CS_MAT_COLS * CS_TWO) as *i64
355 cs_de2000_ex(cs, L1, a1, b1, L2, a2, b2, o)
356 return o[0]
357}
358
359// ---- ICC matrix/TRC profile ----
360func cs_icc_xyz_tag(b: *u8, off: i64, p: *i64, col: i64) -> i64 {
361 let base: i64 = off + CS_ICC_TYPE_HDR_BYTES
362 p[CS_P_MAT + col] = cs_s15f16_q30(b, base)
363 p[CS_P_MAT + CS_MAT_COLS + col] = cs_s15f16_q30(b, base + CS_ICC_XYZ_BYTES)
364 p[CS_P_MAT + CS_MAT_COLS * CS_TWO + col] = cs_s15f16_q30(b, base + CS_ICC_XYZ_BYTES * CS_TWO)
365 return 0
366}
367func cs_icc_wtpt_tag(b: *u8, off: i64, p: *i64) -> i64 {
368 let base: i64 = off + CS_ICC_TYPE_HDR_BYTES
369 p[CS_P_WTPT] = cs_s15f16_q30(b, base)
370 p[CS_P_WTPT + 1] = cs_s15f16_q30(b, base + CS_ICC_XYZ_BYTES)
371 p[CS_P_WTPT + CS_TWO] = cs_s15f16_q30(b, base + CS_ICC_XYZ_BYTES * CS_TWO)
372 return 0
373}
374// curveType -> a 256-entry Q30 linearisation LUT; returns the TRC kind
375func cs_icc_trc_tag(cs: *i64, b: *u8, off: i64, lut: *i64) -> i64 {
376 let fq: *i64 = cs_fq(cs)
377 let count: i64 = cs_be32(b, off + CS_ICC_TYPE_HDR_BYTES)
378 let vals: i64 = off + CS_ICC_TYPE_HDR_BYTES + CS_ICC_XYZ_BYTES
379 var c: i64 = 0
380 if count == 0 {
381 while c < CS_LUT_ENTRIES { lut[c] = fq_div(fq_from_int(c), fq_from_int(CS_BYTE_MAX)); c = c + 1 }
382 return CS_ICC_TRC_IDENTITY
383 }
384 if count == 1 {
385 let gamma: i64 = cs_be16(b, vals) * CS_ICC_U8F8_SCALE
386 while c < CS_LUT_ENTRIES { lut[c] = fq_pow(fq, fq_div(fq_from_int(c), fq_from_int(CS_BYTE_MAX)), gamma); c = c + 1 }
387 return CS_ICC_TRC_GAMMA
388 }
389 while c < CS_LUT_ENTRIES {
390 let num: i64 = c * (count - 1)
391 let idx: i64 = num / CS_BYTE_MAX
392 let rem: i64 = num - idx * CS_BYTE_MAX
393 let v0: i64 = cs_be16(b, vals + idx * CS_TWO)
394 var v1: i64 = v0
395 if idx + 1 < count { v1 = cs_be16(b, vals + (idx + 1) * CS_TWO) }
396 let lin_u16: i64 = v0 * CS_BYTE_MAX + (v1 - v0) * rem
397 // MEASURED 2026-08-24: fq_div(from_int(lin_u16), from_int(65535*255)) handed fq_div a 2^54 denominator, outside
398 // its 2^47 contract; the remainder step overflowed and every table entry above ~0.02 was garbage (white survived
399 // only because a == b). (lin_u16 << 30) / (65535*255) is exact in i64: lin_u16 <= 2^24, so the shift stays < 2^54.
400 lut[c] = (lin_u16 << FQ_SHIFT) / (CS_ICC_U16_MAX * CS_BYTE_MAX)
401 c = c + 1
402 }
403 return CS_ICC_TRC_TABLE
404}
405// 3x3 inverse by adjugate over Q30; returns 0 on a singular matrix
406func cs_inv3(m: *i64, mo: i64, out: *i64, oo: i64) -> i64 {
407 let a: i64 = m[mo]
408 let b: i64 = m[mo + 1]
409 let c: i64 = m[mo + 2]
410 let d: i64 = m[mo + 3]
411 let e: i64 = m[mo + 4]
412 let f: i64 = m[mo + 5]
413 let g: i64 = m[mo + 6]
414 let h: i64 = m[mo + 7]
415 let i: i64 = m[mo + 8]
416 let A: i64 = fq_mul(e, i) - fq_mul(f, h)
417 let B: i64 = fq_mul(f, g) - fq_mul(d, i)
418 let C: i64 = fq_mul(d, h) - fq_mul(e, g)
419 let D: i64 = fq_mul(c, h) - fq_mul(b, i)
420 let E: i64 = fq_mul(a, i) - fq_mul(c, g)
421 let F: i64 = fq_mul(b, g) - fq_mul(a, h)
422 let G: i64 = fq_mul(b, f) - fq_mul(c, e)
423 let H: i64 = fq_mul(c, d) - fq_mul(a, f)
424 let I: i64 = fq_mul(a, e) - fq_mul(b, d)
425 let det: i64 = fq_mul(a, A) + fq_mul(b, B) + fq_mul(c, C)
426 if det == 0 { return 0 }
427 out[oo] = fq_div(A, det)
428 out[oo + 1] = fq_div(D, det)
429 out[oo + 2] = fq_div(G, det)
430 out[oo + 3] = fq_div(B, det)
431 out[oo + 4] = fq_div(E, det)
432 out[oo + 5] = fq_div(H, det)
433 out[oo + 6] = fq_div(C, det)
434 out[oo + 7] = fq_div(F, det)
435 out[oo + 8] = fq_div(I, det)
436 return 1
437}
438// base64 (RFC 4648 alphabet, whitespace skipped, = padding) -> bytes; returns the decoded length. A profile may
439// only reach this estate as TEXT (the MCP write lane carries no binary; the fetchable copy of sRGB2014.icc is
440// Google's git mirror in format=TEXT), so the loader accepts base64 transport when the raw bytes carry no acsp magic.
441func cs_b64_val(c: i64) -> i64 {
442 if c >= CS_B64_UPPER_A { if c <= CS_B64_UPPER_Z { return c - CS_B64_UPPER_A } }
443 if c >= CS_B64_LOWER_A { if c <= CS_B64_LOWER_Z { return c - CS_B64_LOWER_A + CS_B64_LOWER_BASE } }
444 if c >= CS_ASCII_ZERO { if c <= CS_ASCII_NINE { return c - CS_ASCII_ZERO + CS_B64_DIGIT_BASE } }
445 if c == CS_B64_PLUS { return CS_B64_PLUS_VAL }
446 if c == CS_B64_SLASH { return CS_B64_SLASH_VAL }
447 return 0 - 1
448}
449func cs_b64_decode(src: *u8, n: i64, out: *u8) -> i64 {
450 var acc: i64 = 0
451 var bits: i64 = 0
452 var o: i64 = 0
453 var i: i64 = 0
454 while i < n {
455 let c: i64 = src[i] as i64
456 if c == CS_B64_PAD { i = n } else {
457 let v: i64 = cs_b64_val(c)
458 if v >= 0 {
459 acc = (acc << CS_B64_BITS) + v
460 bits = bits + CS_B64_BITS
461 if bits >= CS_BYTE_BITS_1 {
462 bits = bits - CS_BYTE_BITS_1
463 out[o] = ((acc >> bits) & CS_BYTE_MAX) as u8
464 acc = acc & ((1 << bits) - 1)
465 o = o + 1
466 }
467 }
468 i = i + 1
469 }
470 }
471 return o
472}
473// load a matrix/TRC profile (raw ICC bytes, or base64 text of them); p[CS_P_OK] = 1 only when all seven tags parsed
474// and the matrix inverts. p[CS_P_B64] = 1 when the file was base64 transport.
475func cs_icc_load(cs: *i64, path: *u8) -> *i64 {
476 let p: *i64 = sys_mmap(CS_P_SLOTS * CS_SLOT) as *i64
477 p[CS_P_OK] = 0
478 p[CS_P_HAVE] = 0
479 p[CS_P_B64] = 0
480 let lp: *i64 = sys_mmap(CS_SLOT * CS_TWO) as *i64
481 let raw: *u8 = sys_read_file(path, lp)
482 if (raw as i64) == 0 { return p }
483 var b: *u8 = raw
484 var n: i64 = lp[0]
485 if n < CS_ICC_HDR_BYTES + CS_ICC_XYZ_BYTES { return p }
486 if cs_tag_is(b, CS_ICC_MAGIC_OFF, "acsp" as *u8) == 0 {
487 let dec: *u8 = sys_mmap(n)
488 let dn: i64 = cs_b64_decode(raw, n, dec)
489 if dn < CS_ICC_HDR_BYTES + CS_ICC_XYZ_BYTES { p[CS_P_BYTES] = n; return p }
490 if cs_tag_is(dec, CS_ICC_MAGIC_OFF, "acsp" as *u8) == 0 { p[CS_P_BYTES] = n; return p }
491 b = dec
492 n = dn
493 p[CS_P_B64] = 1
494 }
495 p[CS_P_BYTES] = n
496 let ntags: i64 = cs_be32(b, CS_ICC_HDR_BYTES)
497 p[CS_P_TAGS] = ntags
498 let lr: *i64 = sys_mmap(CS_LUT_ENTRIES * CS_SLOT) as *i64
499 let lg: *i64 = sys_mmap(CS_LUT_ENTRIES * CS_SLOT) as *i64
500 let lb: *i64 = sys_mmap(CS_LUT_ENTRIES * CS_SLOT) as *i64
501 p[CS_P_LUT_R] = lr as i64
502 p[CS_P_LUT_G] = lg as i64
503 p[CS_P_LUT_B] = lb as i64
504 var have: i64 = 0
505 var i: i64 = 0
506 while i < ntags {
507 let e: i64 = CS_ICC_HDR_BYTES + CS_ICC_XYZ_BYTES + i * CS_ICC_TAG_ENTRY_BYTES
508 if e + CS_ICC_TAG_ENTRY_BYTES > n { i = ntags } else {
509 let off: i64 = cs_be32(b, e + CS_ICC_TAG_OFF_FIELD)
510 let sz: i64 = cs_be32(b, e + CS_ICC_TAG_SIZE_FIELD)
511 if off + sz <= n { if sz >= CS_ICC_TYPE_HDR_BYTES {
512 if cs_tag_is(b, e, "rXYZ" as *u8) == 1 { cs_icc_xyz_tag(b, off, p, 0); have = have + 1 }
513 if cs_tag_is(b, e, "gXYZ" as *u8) == 1 { cs_icc_xyz_tag(b, off, p, 1); have = have + 1 }
514 if cs_tag_is(b, e, "bXYZ" as *u8) == 1 { cs_icc_xyz_tag(b, off, p, CS_TWO); have = have + 1 }
515 if cs_tag_is(b, e, "wtpt" as *u8) == 1 { cs_icc_wtpt_tag(b, off, p); have = have + 1 }
516 if cs_tag_is(b, e, "rTRC" as *u8) == 1 { p[CS_P_TRC_R] = cs_icc_trc_tag(cs, b, off, lr); have = have + 1 }
517 if cs_tag_is(b, e, "gTRC" as *u8) == 1 { p[CS_P_TRC_G] = cs_icc_trc_tag(cs, b, off, lg); have = have + 1 }
518 if cs_tag_is(b, e, "bTRC" as *u8) == 1 { p[CS_P_TRC_B] = cs_icc_trc_tag(cs, b, off, lb); have = have + 1 }
519 } }
520 i = i + 1
521 }
522 }
523 p[CS_P_HAVE] = have
524 if have == CS_ICC_TAGS_NEEDED { if cs_inv3(p, CS_P_MAT, p, CS_P_INV) == 1 { p[CS_P_OK] = 1 } }
525 // free the FILE buffer only: a base64-decoded copy lives in the arena and must not be munmapped
526 sys_free_file(raw, lp[0])
527 return p
528}
529// 8-bit RGB codes -> Lab (Q30) through the profile (TRC -> matrix -> XYZ D50 -> Lab against wtpt)
530func cs_prof_rgb_to_lab(cs: *i64, p: *i64, r: i64, g: i64, b: i64, out3: *i64) -> i64 {
531 let lr: *i64 = p[CS_P_LUT_R] as *i64
532 let lg: *i64 = p[CS_P_LUT_G] as *i64
533 let lb: *i64 = p[CS_P_LUT_B] as *i64
534 let vr: i64 = lr[r]
535 let vg: i64 = lg[g]
536 let vb: i64 = lb[b]
537 let X: i64 = fq_mul(p[CS_P_MAT], vr) + fq_mul(p[CS_P_MAT + 1], vg) + fq_mul(p[CS_P_MAT + 2], vb)
538 let Y: i64 = fq_mul(p[CS_P_MAT + 3], vr) + fq_mul(p[CS_P_MAT + 4], vg) + fq_mul(p[CS_P_MAT + 5], vb)
539 let Z: i64 = fq_mul(p[CS_P_MAT + 6], vr) + fq_mul(p[CS_P_MAT + 7], vg) + fq_mul(p[CS_P_MAT + 8], vb)
540 return cs_lab_from_xyz(cs, X, Y, Z, p[CS_P_WTPT], p[CS_P_WTPT + 1], p[CS_P_WTPT + 2], out3)
541}
542// inverse LUT: the 8-bit code whose linear value is nearest lin (LUTs are monotone non-decreasing)
543func cs_lut_inverse(lut: *i64, lin: i64) -> i64 {
544 var lo: i64 = 0
545 var hi: i64 = CS_LUT_ENTRIES - 1
546 if lin <= lut[0] { return 0 }
547 if lin >= lut[hi] { return hi }
548 while hi - lo > 1 {
549 let mid: i64 = (lo + hi) / CS_TWO
550 if lut[mid] <= lin { lo = mid } else { hi = mid }
551 }
552 if (lin - lut[lo]) <= (lut[hi] - lin) { return lo }
553 return hi
554}
555func cs_clamp01(v: i64) -> i64 { if v < 0 { return 0 } if v > FQ_ONE { return FQ_ONE } return v }
556// Lab (Q30) -> 8-bit RGB codes through the profile inverse; out3 = r g b codes. Returns 1 if any channel clamped.
557func cs_prof_lab_to_rgb(cs: *i64, p: *i64, L: i64, a: i64, b: i64, out3: *i64) -> i64 {
558 let xyz: *i64 = sys_mmap(CS_SLOT * CS_MAT_COLS) as *i64
559 cs_xyz_from_lab(cs, L, a, b, p[CS_P_WTPT], p[CS_P_WTPT + 1], p[CS_P_WTPT + 2], xyz)
560 let m: i64 = CS_P_INV
561 let r0: i64 = fq_mul(p[m], xyz[0]) + fq_mul(p[m + 1], xyz[1]) + fq_mul(p[m + 2], xyz[2])
562 let g0: i64 = fq_mul(p[m + 3], xyz[0]) + fq_mul(p[m + 4], xyz[1]) + fq_mul(p[m + 5], xyz[2])
563 let b0: i64 = fq_mul(p[m + 6], xyz[0]) + fq_mul(p[m + 7], xyz[1]) + fq_mul(p[m + 8], xyz[2])
564 var clamped: i64 = 0
565 if r0 < 0 { clamped = 1 }
566 if r0 > FQ_ONE { clamped = 1 }
567 if g0 < 0 { clamped = 1 }
568 if g0 > FQ_ONE { clamped = 1 }
569 if b0 < 0 { clamped = 1 }
570 if b0 > FQ_ONE { clamped = 1 }
571 out3[0] = cs_lut_inverse(p[CS_P_LUT_R] as *i64, cs_clamp01(r0))
572 out3[1] = cs_lut_inverse(p[CS_P_LUT_G] as *i64, cs_clamp01(g0))
573 out3[2] = cs_lut_inverse(p[CS_P_LUT_B] as *i64, cs_clamp01(b0))
574 return clamped
575}
576
577// ---- gamut boundary descriptor: segment maxima on a sphere around (L=50, 0, 0) ----
578func cs_gbd_alloc(hs: i64, es: i64) -> *i64 {
579 let g: *i64 = sys_mmap((CS_G_BASE + hs * es) * CS_SLOT) as *i64
580 g[CS_G_H] = hs
581 g[CS_G_E] = es
582 g[CS_G_COUNT] = 0
583 var i: i64 = 0
584 while i < hs * es { g[CS_G_BASE + i] = 0; i = i + 1 }
585 return g
586}
587// out3: [0] radius (Q30), [1] hue sector, [2] elevation sector
588func cs_gbd_sector(cs: *i64, g: *i64, L: i64, a: i64, b: i64, out3: *i64) -> i64 {
589 let fq: *i64 = cs_fq(cs)
590 let dl: i64 = L - fq_from_int(CS_G_CENTER_L)
591 let C: i64 = fq_hypot(a, b)
592 let r: i64 = fq_hypot(dl, C)
593 let h: i64 = cs_hue_deg(fq, a, b)
594 let pi: i64 = fq[FQ_C_PI]
595 let theta: i64 = fq_atan2(fq, dl, C) + pi / CS_TWO
596 var hs: i64 = fq_to_int(fq_div(fq_mul(h, fq_from_int(g[CS_G_H])), fq_from_int(CS_DEG_FULL)))
597 var es: i64 = fq_to_int(fq_div(fq_mul(theta, fq_from_int(g[CS_G_E])), pi))
598 if hs < 0 { hs = 0 }
599 if hs >= g[CS_G_H] { hs = g[CS_G_H] - 1 }
600 if es < 0 { es = 0 }
601 if es >= g[CS_G_E] { es = g[CS_G_E] - 1 }
602 out3[0] = r
603 out3[1] = hs
604 out3[2] = es
605 return 0
606}
607func cs_gbd_add(cs: *i64, g: *i64, L: i64, a: i64, b: i64) -> i64 {
608 let s: *i64 = sys_mmap(CS_SLOT * CS_MAT_COLS) as *i64
609 cs_gbd_sector(cs, g, L, a, b, s)
610 let k: i64 = CS_G_BASE + s[1] * g[CS_G_E] + s[2]
611 if s[0] > g[k] { g[k] = s[0] }
612 g[CS_G_COUNT] = g[CS_G_COUNT] + 1
613 return 0
614}
615// boundary radius of a sector, falling back to the max of its 8 neighbours (hue wraps); CS_G_UNSAMPLED if none
616func cs_gbd_rmax(g: *i64, hs: i64, es: i64) -> i64 {
617 let H: i64 = g[CS_G_H]
618 let E: i64 = g[CS_G_E]
619 let own: i64 = g[CS_G_BASE + hs * E + es]
620 if own > 0 { return own }
621 var best: i64 = 0
622 var dh: i64 = 0 - 1
623 while dh <= 1 {
624 var de: i64 = 0 - 1
625 while de <= 1 {
626 var hh: i64 = hs + dh
627 if hh < 0 { hh = hh + H }
628 if hh >= H { hh = hh - H }
629 let ee: i64 = es + de
630 if ee >= 0 { if ee < E {
631 let v: i64 = g[CS_G_BASE + hh * E + ee]
632 if v > best { best = v }
633 } }
634 de = de + 1
635 }
636 dh = dh + 1
637 }
638 if best > 0 { return best }
639 return CS_G_UNSAMPLED
640}
641func cs_gbd_unsampled(g: *i64) -> i64 {
642 var n: i64 = 0
643 var i: i64 = 0
644 while i < g[CS_G_H] * g[CS_G_E] { if g[CS_G_BASE + i] == 0 { n = n + 1 } i = i + 1 }
645 return n
646}
647// permil of the boundary radius at this point's sector: <= 1000 inside, > 1000 outside, CS_G_UNSAMPLED if the sector
648// and its neighbours were never sampled (an honest third state, never silently inside)
649func cs_gbd_permil(cs: *i64, g: *i64, L: i64, a: i64, b: i64) -> i64 {
650 let s: *i64 = sys_mmap(CS_SLOT * CS_MAT_COLS) as *i64
651 cs_gbd_sector(cs, g, L, a, b, s)
652 let rmax: i64 = cs_gbd_rmax(g, s[1], s[2])
653 if rmax == CS_G_UNSAMPLED { return CS_G_UNSAMPLED }
654 if rmax == 0 { return CS_G_UNSAMPLED }
655 return fq_to_int(fq_div(fq_mul(s[0], fq_from_int(CS_PERMIL)), rmax))
656}
657// sample the six faces of the RGB cube through the profile into the descriptor; returns samples added
658func cs_gbd_from_profile(cs: *i64, p: *i64, g: *i64, steps: i64) -> i64 {
659 let lab: *i64 = sys_mmap(CS_SLOT * CS_MAT_COLS) as *i64
660 var added: i64 = 0
661 var face: i64 = 0
662 while face < CS_CUBE_FACES {
663 var u: i64 = 0
664 while u <= steps {
665 var v: i64 = 0
666 while v <= steps {
667 let cu: i64 = (u * CS_BYTE_MAX) / steps
668 let cv: i64 = (v * CS_BYTE_MAX) / steps
669 var r: i64 = 0
670 var gg: i64 = 0
671 var bb: i64 = 0
672 if face == 0 { r = 0; gg = cu; bb = cv }
673 if face == 1 { r = CS_BYTE_MAX; gg = cu; bb = cv }
674 if face == 2 { r = cu; gg = 0; bb = cv }
675 if face == 3 { r = cu; gg = CS_BYTE_MAX; bb = cv }
676 if face == 4 { r = cu; gg = cv; bb = 0 }
677 if face == 5 { r = cu; gg = cv; bb = CS_BYTE_MAX }
678 cs_prof_rgb_to_lab(cs, p, r, gg, bb, lab)
679 cs_gbd_add(cs, g, lab[0], lab[1], lab[2])
680 added = added + 1
681 v = v + 1
682 }
683 u = u + 1
684 }
685 face = face + 1
686 }
687 return added
688}
689// cs_icc_gamut: THE SUBSTRATE GAMUT FROM A PROFILE in one call (the PT1 watch contract): load the matrix/TRC profile
690// and sample its cube faces into a fresh descriptor. Returns 0 when the profile does not parse (the caller refuses).
691func cs_icc_gamut(cs: *i64, path: *u8, hs: i64, es: i64, steps: i64) -> *i64 {
692 let p: *i64 = cs_icc_load(cs, path)
693 if p[CS_P_OK] != 1 { return 0 as *i64 }
694 let g: *i64 = cs_gbd_alloc(hs, es)
695 cs_gbd_from_profile(cs, p, g, steps)
696 return g
697}
698// CHROMA CLIPPING at constant L and hue: bisect the chroma scale t in [0,1] to the boundary. out3 = mapped Lab.
699// Returns the pre-clip permil (<= 1000 means the point was already inside and is returned unchanged).
700func cs_gamut_clip(cs: *i64, g: *i64, L: i64, a: i64, b: i64, iters: i64, out3: *i64) -> i64 {
701 let before: i64 = cs_gbd_permil(cs, g, L, a, b)
702 out3[0] = L
703 out3[1] = a
704 out3[2] = b
705 if before == CS_G_UNSAMPLED { return before }
706 if before <= CS_PERMIL { return before }
707 var lo: i64 = CS_CLIP_LOW_RADIUS
708 var hi: i64 = FQ_ONE
709 var i: i64 = 0
710 while i < iters {
711 let mid: i64 = (lo + hi) / CS_TWO
712 let pm: i64 = cs_gbd_permil(cs, g, L, fq_mul(a, mid), fq_mul(b, mid))
713 var inside: i64 = 0
714 if pm != CS_G_UNSAMPLED { if pm <= CS_PERMIL { inside = 1 } }
715 if inside == 1 { lo = mid } else { hi = mid }
716 i = i + 1
717 }
718 out3[1] = fq_mul(a, lo)
719 out3[2] = fq_mul(b, lo)
720 return before
721}
722
723// ---- EXACT GAMUT TEST FOR A MATRIX/TRC PROFILE (2026-08-24) ----
724// A colour is inside a matrix/TRC gamut iff its LINEAR RGB through the profile inverse lies in [0,1] on every channel:
725// the cube IS the gamut, no descriptor approximation is needed. MEASURED the day this shipped: the 16-step segment-
726// maxima descriptor left 536 of 648 sphere sectors unsampled and read 18 percent of an sRGB photograph OUTSIDE sRGB;
727// this test reads 0 by construction. The descriptor stays for MEASURED substrates (swatch Lab samples have no matrix).
728// Lab -> linear RGB (Q30) through the inverse matrix, unclamped; out3 = r g b linear
729func cs_prof_lab_to_lin(cs: *i64, p: *i64, L: i64, a: i64, b: i64, out3: *i64) -> i64 {
730 let xyz: *i64 = sys_mmap(CS_SLOT * CS_MAT_COLS) as *i64
731 cs_xyz_from_lab(cs, L, a, b, p[CS_P_WTPT], p[CS_P_WTPT + 1], p[CS_P_WTPT + 2], xyz)
732 let m: i64 = CS_P_INV
733 out3[0] = fq_mul(p[m], xyz[0]) + fq_mul(p[m + 1], xyz[1]) + fq_mul(p[m + 2], xyz[2])
734 out3[1] = fq_mul(p[m + 3], xyz[0]) + fq_mul(p[m + 4], xyz[1]) + fq_mul(p[m + 5], xyz[2])
735 out3[2] = fq_mul(p[m + 6], xyz[0]) + fq_mul(p[m + 7], xyz[1]) + fq_mul(p[m + 8], xyz[2])
736 return 0
737}
738// permil of the linear range: 1000 on the boundary, below 1000 inside (1000 minus the tightest headroom), above 1000
739// outside (1000 plus the largest excursion). Never UNSAMPLED: a matrix gamut has no holes.
740func cs_prof_gamut_permil(cs: *i64, p: *i64, L: i64, a: i64, b: i64) -> i64 {
741 let lin: *i64 = sys_mmap(CS_SLOT * CS_MAT_COLS) as *i64
742 cs_prof_lab_to_lin(cs, p, L, a, b, lin)
743 var excess: i64 = 0 - FQ_ONE
744 var i: i64 = 0
745 while i < CS_MAT_COLS {
746 let lo: i64 = 0 - lin[i]
747 let hi: i64 = lin[i] - FQ_ONE
748 if lo > excess { excess = lo }
749 if hi > excess { excess = hi }
750 i = i + 1
751 }
752 return CS_PERMIL + fq_to_int(fq_mul(excess, fq_from_int(CS_PERMIL)))
753}
754// ---- CHROMATIC ADAPTATION ----
755// one Q30 coefficient from a published matrix carried at 1/10000
756func cs_cat_q(v: i64) -> i64 { return fq_div(fq_from_int(v), fq_from_int(CS_CAT_SCALE)) }
757// the cone-response matrix for the chosen transform, Q30, row-major in out9. Returns the transform id it
758// actually wrote, or -1 for an unknown id: an unrecognised transform REFUSES rather than defaulting, because
759// a silent fallback to Bradford would make a CAT02 request that never happened read as though it had.
760func cs_cat_cone(which: i64, out9: *i64) -> i64 {
761 if which == CS_CAT_BRADFORD {
762 out9[0] = cs_cat_q(CS_CAT_BFD_0)
763 out9[1] = cs_cat_q(CS_CAT_BFD_1)
764 out9[2] = cs_cat_q(CS_CAT_BFD_2)
765 out9[3] = cs_cat_q(CS_CAT_BFD_3)
766 out9[4] = cs_cat_q(CS_CAT_BFD_4)
767 out9[5] = cs_cat_q(CS_CAT_BFD_5)
768 out9[6] = cs_cat_q(CS_CAT_BFD_6)
769 out9[7] = cs_cat_q(CS_CAT_BFD_7)
770 out9[8] = cs_cat_q(CS_CAT_BFD_8)
771 return CS_CAT_BRADFORD
772 }
773 if which == CS_CAT_CAT02 {
774 out9[0] = cs_cat_q(CS_CAT_02_0)
775 out9[1] = cs_cat_q(CS_CAT_02_1)
776 out9[2] = cs_cat_q(CS_CAT_02_2)
777 out9[3] = cs_cat_q(CS_CAT_02_3)
778 out9[4] = cs_cat_q(CS_CAT_02_4)
779 out9[5] = cs_cat_q(CS_CAT_02_5)
780 out9[6] = cs_cat_q(CS_CAT_02_6)
781 out9[7] = cs_cat_q(CS_CAT_02_7)
782 out9[8] = cs_cat_q(CS_CAT_02_8)
783 return CS_CAT_CAT02
784 }
785 return 0 - 1
786}
787// out9 = a * b, 3x3 Q30
788func cs_cat_mul3(a: *i64, b: *i64, out9: *i64) -> i64 {
789 var r: i64 = 0
790 while r < CS_MAT_COLS {
791 var c: i64 = 0
792 while c < CS_MAT_COLS {
793 var s: i64 = 0
794 var k: i64 = 0
795 while k < CS_MAT_COLS {
796 s = s + fq_mul(a[r * CS_MAT_COLS + k], b[k * CS_MAT_COLS + c])
797 k = k + 1
798 }
799 out9[r * CS_MAT_COLS + c] = s
800 c = c + 1
801 }
802 r = r + 1
803 }
804 return 0
805}
806// out3 = m9 * (X,Y,Z)
807func cs_cat_apply(m9: *i64, X: i64, Y: i64, Z: i64, out3: *i64) -> i64 {
808 out3[0] = fq_mul(m9[0], X) + fq_mul(m9[1], Y) + fq_mul(m9[2], Z)
809 out3[1] = fq_mul(m9[3], X) + fq_mul(m9[4], Y) + fq_mul(m9[5], Z)
810 out3[2] = fq_mul(m9[6], X) + fq_mul(m9[7], Y) + fq_mul(m9[8], Z)
811 return 0
812}
813// the full adaptation matrix Minv * D * M, taking XYZ measured under the SOURCE white to XYZ under the
814// DESTINATION white. Returns 1 on success, 0 on refusal (unknown transform id, singular cone matrix, or a
815// zero source cone response). A refusal NEVER writes out9, because an unadapted matrix that reads adapted
816// is precisely the defect this function exists to remove.
817func cs_cat_build(Xs: i64, Ys: i64, Zs: i64, Xd: i64, Yd: i64, Zd: i64, which: i64, out9: *i64) -> i64 {
818 let M: *i64 = sys_mmap(CS_CAT_M * CS_SLOT) as *i64
819 let Mi: *i64 = sys_mmap(CS_CAT_M * CS_SLOT) as *i64
820 let D: *i64 = sys_mmap(CS_CAT_M * CS_SLOT) as *i64
821 let T: *i64 = sys_mmap(CS_CAT_M * CS_SLOT) as *i64
822 let ws: *i64 = sys_mmap(CS_MAT_COLS * CS_SLOT) as *i64
823 let wd: *i64 = sys_mmap(CS_MAT_COLS * CS_SLOT) as *i64
824 if cs_cat_cone(which, M) < 0 { return 0 }
825 if cs_inv3(M, 0, Mi, 0) != 1 { return 0 }
826 cs_cat_apply(M, Xs, Ys, Zs, ws)
827 cs_cat_apply(M, Xd, Yd, Zd, wd)
828 if ws[0] == 0 { return 0 }
829 if ws[1] == 0 { return 0 }
830 if ws[2] == 0 { return 0 }
831 var i: i64 = 0
832 while i < CS_CAT_M {
833 D[i] = 0
834 i = i + 1
835 }
836 D[0] = fq_div(wd[0], ws[0])
837 D[4] = fq_div(wd[1], ws[1])
838 D[8] = fq_div(wd[2], ws[2])
839 cs_cat_mul3(D, M, T)
840 cs_cat_mul3(Mi, T, out9)
841 return 1
842}
843// adapt one XYZ from the source white to the destination white in one call. Returns 1 on success, 0 on a
844// refusal from cs_cat_build, and on refusal out3 is left UNTOUCHED so a caller cannot read a failed
845// adaptation as an identity one.
846func cs_cat_xyz(X: i64, Y: i64, Z: i64, Xs: i64, Ys: i64, Zs: i64, Xd: i64, Yd: i64, Zd: i64, which: i64, out3: *i64) -> i64 {
847 let m: *i64 = sys_mmap(CS_CAT_M * CS_SLOT) as *i64
848 if cs_cat_build(Xs, Ys, Zs, Xd, Yd, Zd, which, m) != 1 { return 0 }
849 cs_cat_apply(m, X, Y, Z, out3)
850 return 1
851}
852
853// chroma clipping at constant L and hue against the EXACT matrix gamut; returns the pre-clip permil, out3 = mapped Lab
854func cs_prof_gamut_clip(cs: *i64, p: *i64, L: i64, a: i64, b: i64, iters: i64, out3: *i64) -> i64 {
855 let before: i64 = cs_prof_gamut_permil(cs, p, L, a, b)
856 out3[0] = L
857 out3[1] = a
858 out3[2] = b
859 if before <= CS_PERMIL { return before }
860 var lo: i64 = CS_CLIP_LOW_RADIUS
861 var hi: i64 = FQ_ONE
862 var i: i64 = 0
863 while i < iters {
864 let mid: i64 = (lo + hi) / CS_TWO
865 let pm: i64 = cs_prof_gamut_permil(cs, p, L, fq_mul(a, mid), fq_mul(b, mid))
866 if pm <= CS_PERMIL { lo = mid } else { hi = mid }
867 i = i + 1
868 }
869 out3[1] = fq_mul(a, lo)
870 out3[2] = fq_mul(b, lo)
871 return before
872}