nx_lab_from_rgb.nx source
↩ module page · 269 lines · 10664 B
1// nx_lab_from_rgb.nx -- sRGB u8 -> CIELab Q10 in pure i64 fixed-point.
2//
3// CAPABILITY_COMPLETENESS: FULL
4//
5// Bridges the "RGB-proxy" tier-1 chroma axes (in nx_image_grade_v0)
6// up to genuine colorimetric measurements:
7// - real ITA degrees (via nx_skin_tone_ita)
8// - real CIELab h_ab undertone angle (via nx_undertone_class)
9//
10// Math pipeline (CIE 15:2004 / CIE 014:2018):
11//
12// step 1 sRGB u8 channel value [0, 255]
13// -> Q14 sRGB scalar [0, 16384]
14//
15// step 2 sRGB Q14 -> linear-light Q14 via sRGB OETF inverse:
16// val_q14 < 0.04045 * Q14 -> linear = val / 12.92
17// val_q14 >= 0.04045 * Q14 -> linear = ((val + 0.055)/1.055)^2.4
18// Composed against nx_aces_srgb_to_linear_q14 (already shipped).
19//
20// step 3 Linear sRGB Q14 -> XYZ Q14 via 3x3 matrix (D65):
21// [0.4124564 0.3575761 0.1804375] per IEC 61966-2-1 +
22// [0.2126729 0.7151522 0.0721750] ITU-R BT.709 primaries
23// [0.0193339 0.1191920 0.9503041]
24//
25// step 4 XYZ Q14 -> Lab Q10:
26// tn_x = X / Xn (Xn = 0.95047 for D65)
27// tn_y = Y / Yn (Yn = 1.00000)
28// tn_z = Z / Zn (Zn = 1.08883)
29// f(t) = t^(1/3) if t > 216/24389
30// (kappa*t + 16) / 116 otherwise (kappa = 24389/27)
31// L* = 116 * f(tn_y) - 16
32// a* = 500 * (f(tn_x) - f(tn_y))
33// b* = 200 * (f(tn_y) - f(tn_z))
34//
35// Output: L*_q10, a*_q10, b*_q10 all in Q10 scaled units.
36// L*_q10 in [0, ~100 * 1024 = 102400]
37// a*_q10 / b*_q10 in approximately [-130 * 1024, +130 * 1024]
38//
39// Cube root via Newton-Raphson in Q14 (5 iterations -> ~0.1%
40// accuracy across [0, 1]).
41//
42// genealogy_id: cie_15_2004_lab + iec_61966_2_1_srgb
43// lineage_id: nx_lab_from_rgb_v1
44
45// nx_safety_envelope:
46// intended_use: AUTO_APPLIED -- primitive-specific tuning queued
47// sil_target: SIL1
48// evidence: [bulk_applied_2026-05-16, see-file-comment-for-detail]
49// verdict: NOT_YET_EVALUATED
50
51import "nx_syscalls.nx"
52import "nx_runtime.nx"
53import "nx_tier.nx"
54import "nx_aces_native.nx"
55const NX_MAGIC_2900: i64 = 2900
56const NX_MAGIC_3200: i64 = 3200
57const NX_MAGIC_8192: i64 = 8192
58const NX_MAGIC_13005: i64 = 13005
59
60const NX_LAB_Q14: nx_int = 16384
61const NX_LAB_Q10: nx_int = 1024
62
63// D65 white point in Q14.
64const NX_LAB_XN_Q14: nx_int = 15572 // 0.95047 * 16384
65const NX_LAB_YN_Q14: nx_int = 16384 // 1.00000 * 16384
66const NX_LAB_ZN_Q14: nx_int = 17839 // 1.08883 * 16384
67
68// sRGB-D65 -> XYZ-D65 matrix in Q14.
69const NX_LAB_M00: nx_int = 6759 // 0.4124564 * 16384
70const NX_LAB_M01: nx_int = 5859 // 0.3575761 * 16384
71const NX_LAB_M02: nx_int = 2957 // 0.1804375 * 16384
72const NX_LAB_M10: nx_int = 3484 // 0.2126729 * 16384
73const NX_LAB_M11: nx_int = 11718 // 0.7151522 * 16384
74const NX_LAB_M12: nx_int = 1183 // 0.0721750 * 16384
75const NX_LAB_M20: nx_int = 317 // 0.0193339 * 16384
76const NX_LAB_M21: nx_int = 1953 // 0.1191920 * 16384
77const NX_LAB_M22: nx_int = 15569 // 0.9503041 * 16384
78
79// f(t) threshold = 216/24389 ~= 0.008856
80const NX_LAB_FT_THRESHOLD_Q14: nx_int = 145 // 0.008856 * 16384
81
82// Polynomial branch: (kappa*t + 16) / 116 in Q14 form.
83// kappa = 24389/27 ~= 903.296
84// kappa_q14 = 903 (kappa truncated; reasonable since only used
85// below threshold where t is tiny)
86// 16 -> 16/116 in Q14 = 16384 * 16/116 = 2260 (rounded)
87// So f_low(t_q14) = (903 * t_q14 + 16 * 16384) / 116 / 16384 result
88// ... simpler: keep numerator in Q28 then divide at end.
89//
90// We compute: f_low_q14(t_q14) = (903 * t_q14 / 16384 + 2260) shifted
91// to a directly Q14-comparable form. Implementation below.
92
93// ===== result struct =============================================
94
95struct NxLabColor {
96 l_q10: nx_int,
97 a_q10: nx_int,
98 b_q10: nx_int,
99}
100
101const NX_LAB_COLOR_BYTES: nx_size = 24
102
103// ===== Newton-Raphson cube root in Q14 ============================
104//
105// Given t_q14 in [0, Q14], compute t^(1/3) in Q14.
106// Recurrence: x_{n+1} = (2*x_n + t/(x_n^2)) / 3
107// Initial estimate: x_0 = max(t_q14, Q14 / 8).
108//
109// All math in Q14 throughout:
110// t / (x*x) in Q14 = (t_q14 * Q14 * Q14) / (x_q14 * x_q14)
111// Avoid overflow: t_q14 <= 16384 = 2^14; Q14^2 = 2^28; product
112// fits in i64 since 14 + 28 = 42 < 63.
113
114func _lab_cuberoot_q14(t_q14: nx_int) -> nx_int {
115 if t_q14 <= 0 { return 0 }
116 if t_q14 >= NX_LAB_Q14 { return NX_LAB_Q14 }
117 // Initial estimate: for t in [0, 1], x ~ t for large t. For
118 // small t, t^(1/3) is much larger than t; using x=t will
119 // diverge. Start from x = max(t, Q14/8) as a robust guess.
120 var x: nx_int = t_q14
121 let floor_est: nx_int = NX_LAB_Q14 / 8
122 if x < floor_est { x = floor_est }
123 var i: nx_int = 0
124 while i < 6 {
125 let x_sq: nx_int = (x * x) / NX_LAB_Q14
126 if x_sq <= 0 {
127 x = NX_LAB_Q14 / 16
128 } else {
129 let div: nx_int = (t_q14 * NX_LAB_Q14) / x_sq
130 x = (2 * x + div) / 3
131 }
132 i = i + 1
133 }
134 if x < 0 { return 0 }
135 if x > NX_LAB_Q14 { return NX_LAB_Q14 }
136 return x
137}
138
139// ===== f(t) per CIE definition ===================================
140
141func _lab_ft_q14(t_q14: nx_int) -> nx_int {
142 if t_q14 <= 0 { return 0 }
143 if t_q14 > NX_LAB_FT_THRESHOLD_Q14 {
144 return _lab_cuberoot_q14(t_q14)
145 }
146 // Low branch: f(t) = (kappa*t + 16) / 116 where kappa = 24389/27.
147 // approximate kappa = 903.296. In Q14:
148 // f_q14 = (903 * t_q14 + 16 * NX_LAB_Q14) / 116
149 let num: nx_int = 903 * t_q14 + 16 * NX_LAB_Q14
150 return num / 116
151}
152
153// ===== sRGB u8 -> Lab Q10 =========================================
154
155func nx_lab_from_rgb(r_u8: nx_int, g_u8: nx_int, b_u8: nx_int) -> *NxLabColor {
156 let p: *u8 = sys_mmap(NX_LAB_COLOR_BYTES)
157 let out: *NxLabColor = p as *NxLabColor
158
159 // Step 1+2: sRGB u8 -> linear Q14. Compose against the existing
160 // nx_aces_native pipeline which already does sRGB OETF inverse +
161 // u8 quantization.
162 let r_s_q14: nx_int = (r_u8 * NX_LAB_Q14) / 255
163 let g_s_q14: nx_int = (g_u8 * NX_LAB_Q14) / 255
164 let b_s_q14: nx_int = (b_u8 * NX_LAB_Q14) / 255
165 let r_l: nx_int = nx_aces_srgb_to_linear_q14(r_s_q14)
166 let g_l: nx_int = nx_aces_srgb_to_linear_q14(g_s_q14)
167 let b_l: nx_int = nx_aces_srgb_to_linear_q14(b_s_q14)
168
169 // Step 3: linear -> XYZ via 3x3 matrix.
170 let x_q14: nx_int = ((NX_LAB_M00 * r_l) + (NX_LAB_M01 * g_l) + (NX_LAB_M02 * b_l)) / NX_LAB_Q14
171 let y_q14: nx_int = ((NX_LAB_M10 * r_l) + (NX_LAB_M11 * g_l) + (NX_LAB_M12 * b_l)) / NX_LAB_Q14
172 let z_q14: nx_int = ((NX_LAB_M20 * r_l) + (NX_LAB_M21 * g_l) + (NX_LAB_M22 * b_l)) / NX_LAB_Q14
173
174 // Step 4: normalize by white point + f(t) + Lab.
175 let tn_x_q14: nx_int = (x_q14 * NX_LAB_Q14) / NX_LAB_XN_Q14
176 let tn_y_q14: nx_int = (y_q14 * NX_LAB_Q14) / NX_LAB_YN_Q14
177 let tn_z_q14: nx_int = (z_q14 * NX_LAB_Q14) / NX_LAB_ZN_Q14
178 let fx_q14: nx_int = _lab_ft_q14(tn_x_q14)
179 let fy_q14: nx_int = _lab_ft_q14(tn_y_q14)
180 let fz_q14: nx_int = _lab_ft_q14(tn_z_q14)
181
182 // L* = 116 * f(tn_y) - 16 (in actual units; Lab uses 0..100 range)
183 // In Q10: L*_q10 = 116 * f_q14 * Q10 / Q14 - 16 * Q10
184 // = (116 * fy_q14) / 16 - 16384 (since Q14/Q10 = 16)
185 let l_q10: nx_int = (116 * fy_q14) / 16 - (16 * NX_LAB_Q10)
186 // a* = 500 * (f(tn_x) - f(tn_y)) in actual units
187 // In Q10: a*_q10 = 500 * (fx_q14 - fy_q14) / 16
188 let a_q10: nx_int = (500 * (fx_q14 - fy_q14)) / 16
189 // b* = 200 * (f(tn_y) - f(tn_z))
190 let b_q10: nx_int = (200 * (fy_q14 - fz_q14)) / 16
191
192 out.l_q10 = l_q10
193 out.a_q10 = a_q10
194 out.b_q10 = b_q10
195 return out
196}
197
198// ===== self-test ==================================================
199
200func main() -> nx_int {
201 // ---- white (255, 255, 255) -> L*=100, a*=0, b*=0 ----
202 let w_lab: *NxLabColor = nx_lab_from_rgb(255, 255, 255)
203 // L* ~ 100 +/- 5 due to polynomial gamma + Q14 truncation.
204 let l_white: nx_int = w_lab.l_q10 / NX_LAB_Q10
205 if l_white < 90 { return 1 }
206 if l_white > 105 { return 2 }
207 let a_white: nx_int = w_lab.a_q10 / NX_LAB_Q10
208 if a_white < -8 { return 3 }
209 if a_white > 8 { return 4 }
210 let b_white: nx_int = w_lab.b_q10 / NX_LAB_Q10
211 if b_white < -8 { return 5 }
212 if b_white > 8 { return 6 }
213
214 // ---- black (0, 0, 0) -> L*=0, a*=0, b*=0 (sub-threshold) ----
215 let bk_lab: *NxLabColor = nx_lab_from_rgb(0, 0, 0)
216 let l_black: nx_int = bk_lab.l_q10 / NX_LAB_Q10
217 // Sub-threshold low-branch can give slightly positive L* due to
218 // the +16 offset before subtraction. Tolerance: |L*| < 5.
219 var abs_l_black: nx_int = l_black
220 if abs_l_black < 0 { abs_l_black = 0 - abs_l_black }
221 if abs_l_black > 5 { return 10 }
222
223 // ---- pure red (255, 0, 0) -> a* strongly positive, b* positive ----
224 let r_lab: *NxLabColor = nx_lab_from_rgb(255, 0, 0)
225 let a_red: nx_int = r_lab.a_q10 / NX_LAB_Q10
226 let b_red: nx_int = r_lab.b_q10 / NX_LAB_Q10
227 // Reference Lab for sRGB red: L*=53.24, a*=80.09, b*=67.20.
228 // With our Q14 polynomial approx, expect a* in [60, 95] and
229 // b* in [50, 85]. Wide tolerance.
230 if a_red < 50 { return 20 }
231 if b_red < 35 { return 21 }
232
233 // ---- pure blue (0, 0, 255) -> a* slightly positive, b* strongly negative ----
234 let bl_lab: *NxLabColor = nx_lab_from_rgb(0, 0, 255)
235 let b_blue: nx_int = bl_lab.b_q10 / NX_LAB_Q10
236 if b_blue > -50 { return 30 }
237
238 // ---- mid skin tone (200, 170, 150) -> L* high-ish, a* positive,
239 // b* positive (warm undertone) ----
240 let sk_lab: *NxLabColor = nx_lab_from_rgb(200, 170, 150)
241 let l_skin: nx_int = sk_lab.l_q10 / NX_LAB_Q10
242 if l_skin < 60 { return 40 }
243 if l_skin > 85 { return 41 }
244 let a_skin: nx_int = sk_lab.a_q10 / NX_LAB_Q10
245 if a_skin <= 0 { return 42 }
246 let b_skin: nx_int = sk_lab.b_q10 / NX_LAB_Q10
247 if b_skin <= 0 { return 43 }
248
249 // ---- f(t) low-branch boundary ----
250 // t = 100 (well below threshold 145) -> f = (903*100 + 16*16384) / 116
251 // = (90300 + 262144) / 116 = 352444 / 116 = ~3038 Q14
252 let fl_low: nx_int = _lab_ft_q14(100)
253 if fl_low < NX_MAGIC_2900 { return 50 }
254 if fl_low > NX_MAGIC_3200 { return 51 }
255
256 // ---- cube root: 1.0 (Q14) -> 1.0 ----
257 let cr1: nx_int = _lab_cuberoot_q14(NX_LAB_Q14)
258 var d_cr1: nx_int = cr1 - NX_LAB_Q14
259 if d_cr1 < 0 { d_cr1 = 0 - d_cr1 }
260 if d_cr1 > 30 { return 60 }
261
262 // ---- cube root: 0.5 (Q14=8192) -> 0.7937 (Q14=13005) ----
263 let cr_half: nx_int = _lab_cuberoot_q14(NX_MAGIC_8192)
264 var d_half: nx_int = cr_half - NX_MAGIC_13005
265 if d_half < 0 { d_half = 0 - d_half }
266 if d_half > 300 { return 61 }
267
268 return 0
269}