1// Copyright(c) 2021 Björn Ottosson
2//
3// Permission is hereby granted, free of charge, to any person obtaining a copy of
4// this software and associated documentation files(the "Software"), to deal in
5// the Software without restriction, including without limitation the rights to
6// use, copy, modify, merge, publish, distribute, sublicense, and /or sell copies
7// of the Software, and to permit persons to whom the Software is furnished to do
8// so, subject to the following conditions :
9// The above copyright notice and this permission notice shall be included in all
10// copies or substantial portions of the Software.
11// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
12// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
13// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE
14// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
15// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
16// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
17// SOFTWARE.
18
19#ifndef OK_COLOR_H
20#define OK_COLOR_H
21
22#include <cmath>
23#include <cfloat>
24
25class ok_color
26{
27public:
28
29struct Lab { float L; float a; float b; };
30struct RGB { float r; float g; float b; };
31struct HSV { float h; float s; float v; };
32struct HSL { float h; float s; float l; };
33struct LC { float L; float C; };
34
35// Alternative representation of (L_cusp, C_cusp)
36// Encoded so S = C_cusp/L_cusp and T = C_cusp/(1-L_cusp)
37// The maximum value for C in the triangle is then found as fmin(S*L, T*(1-L)), for a given L
38struct ST { float S; float T; };
39
40static constexpr float pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062f;
41
42static float clamp(float x, float min, float max)
43{
44 if (x < min)
45 return min;
46 if (x > max)
47 return max;
48
49 return x;
50}
51
52static float sgn(float x)
53{
54 return (float)(0.f < x) - (float)(x < 0.f);
55}
56
57static float srgb_transfer_function(float a)
58{
59 return .0031308f >= a ? 12.92f * a : 1.055f * powf(a, .4166666666666667f) - .055f;
60}
61
62static float srgb_transfer_function_inv(float a)
63{
64 return .04045f < a ? powf((a + .055f) / 1.055f, 2.4f) : a / 12.92f;
65}
66
67static Lab linear_srgb_to_oklab(RGB c)
68{
69 float l = 0.4122214708f * c.r + 0.5363325363f * c.g + 0.0514459929f * c.b;
70 float m = 0.2119034982f * c.r + 0.6806995451f * c.g + 0.1073969566f * c.b;
71 float s = 0.0883024619f * c.r + 0.2817188376f * c.g + 0.6299787005f * c.b;
72
73 float l_ = cbrtf(l);
74 float m_ = cbrtf(m);
75 float s_ = cbrtf(s);
76
77 return {
78 0.2104542553f * l_ + 0.7936177850f * m_ - 0.0040720468f * s_,
79 1.9779984951f * l_ - 2.4285922050f * m_ + 0.4505937099f * s_,
80 0.0259040371f * l_ + 0.7827717662f * m_ - 0.8086757660f * s_,
81 };
82}
83
84static RGB oklab_to_linear_srgb(Lab c)
85{
86 float l_ = c.L + 0.3963377774f * c.a + 0.2158037573f * c.b;
87 float m_ = c.L - 0.1055613458f * c.a - 0.0638541728f * c.b;
88 float s_ = c.L - 0.0894841775f * c.a - 1.2914855480f * c.b;
89
90 float l = l_ * l_ * l_;
91 float m = m_ * m_ * m_;
92 float s = s_ * s_ * s_;
93
94 return {
95 +4.0767416621f * l - 3.3077115913f * m + 0.2309699292f * s,
96 -1.2684380046f * l + 2.6097574011f * m - 0.3413193965f * s,
97 -0.0041960863f * l - 0.7034186147f * m + 1.7076147010f * s,
98 };
99}
100
101// Finds the maximum saturation possible for a given hue that fits in sRGB
102// Saturation here is defined as S = C/L
103// a and b must be normalized so a^2 + b^2 == 1
104static float compute_max_saturation(float a, float b)
105{
106 // Max saturation will be when one of r, g or b goes below zero.
107
108 // Select different coefficients depending on which component goes below zero first
109 float k0, k1, k2, k3, k4, wl, wm, ws;
110
111 if (-1.88170328f * a - 0.80936493f * b > 1)
112 {
113 // Red component
114 k0 = +1.19086277f; k1 = +1.76576728f; k2 = +0.59662641f; k3 = +0.75515197f; k4 = +0.56771245f;
115 wl = +4.0767416621f; wm = -3.3077115913f; ws = +0.2309699292f;
116 }
117 else if (1.81444104f * a - 1.19445276f * b > 1)
118 {
119 // Green component
120 k0 = +0.73956515f; k1 = -0.45954404f; k2 = +0.08285427f; k3 = +0.12541070f; k4 = +0.14503204f;
121 wl = -1.2684380046f; wm = +2.6097574011f; ws = -0.3413193965f;
122 }
123 else
124 {
125 // Blue component
126 k0 = +1.35733652f; k1 = -0.00915799f; k2 = -1.15130210f; k3 = -0.50559606f; k4 = +0.00692167f;
127 wl = -0.0041960863f; wm = -0.7034186147f; ws = +1.7076147010f;
128 }
129
130 // Approximate max saturation using a polynomial:
131 float S = k0 + k1 * a + k2 * b + k3 * a * a + k4 * a * b;
132
133 // Do one step Halley's method to get closer
134 // this gives an error less than 10e6, except for some blue hues where the dS/dh is close to infinite
135 // this should be sufficient for most applications, otherwise do two/three steps
136
137 float k_l = +0.3963377774f * a + 0.2158037573f * b;
138 float k_m = -0.1055613458f * a - 0.0638541728f * b;
139 float k_s = -0.0894841775f * a - 1.2914855480f * b;
140
141 {
142 float l_ = 1.f + S * k_l;
143 float m_ = 1.f + S * k_m;
144 float s_ = 1.f + S * k_s;
145
146 float l = l_ * l_ * l_;
147 float m = m_ * m_ * m_;
148 float s = s_ * s_ * s_;
149
150 float l_dS = 3.f * k_l * l_ * l_;
151 float m_dS = 3.f * k_m * m_ * m_;
152 float s_dS = 3.f * k_s * s_ * s_;
153
154 float l_dS2 = 6.f * k_l * k_l * l_;
155 float m_dS2 = 6.f * k_m * k_m * m_;
156 float s_dS2 = 6.f * k_s * k_s * s_;
157
158 float f = wl * l + wm * m + ws * s;
159 float f1 = wl * l_dS + wm * m_dS + ws * s_dS;
160 float f2 = wl * l_dS2 + wm * m_dS2 + ws * s_dS2;
161
162 S = S - f * f1 / (f1 * f1 - 0.5f * f * f2);
163 }
164
165 return S;
166}
167
168// finds L_cusp and C_cusp for a given hue
169// a and b must be normalized so a^2 + b^2 == 1
170static LC find_cusp(float a, float b)
171{
172 // First, find the maximum saturation (saturation S = C/L)
173 float S_cusp = compute_max_saturation(a, b);
174
175 // Convert to linear sRGB to find the first point where at least one of r,g or b >= 1:
176 RGB rgb_at_max = oklab_to_linear_srgb({ 1, S_cusp * a, S_cusp * b });
177 float L_cusp = cbrtf(1.f / fmax(fmax(rgb_at_max.r, rgb_at_max.g), rgb_at_max.b));
178 float C_cusp = L_cusp * S_cusp;
179
180 return { L_cusp , C_cusp };
181}
182
183// Finds intersection of the line defined by
184// L = L0 * (1 - t) + t * L1;
185// C = t * C1;
186// a and b must be normalized so a^2 + b^2 == 1
187static float find_gamut_intersection(float a, float b, float L1, float C1, float L0, LC cusp)
188{
189 // Find the intersection for upper and lower half seprately
190 float t;
191 if (((L1 - L0) * cusp.C - (cusp.L - L0) * C1) <= 0.f)
192 {
193 // Lower half
194
195 t = cusp.C * L0 / (C1 * cusp.L + cusp.C * (L0 - L1));
196 }
197 else
198 {
199 // Upper half
200
201 // First intersect with triangle
202 t = cusp.C * (L0 - 1.f) / (C1 * (cusp.L - 1.f) + cusp.C * (L0 - L1));
203
204 // Then one step Halley's method
205 {
206 float dL = L1 - L0;
207 float dC = C1;
208
209 float k_l = +0.3963377774f * a + 0.2158037573f * b;
210 float k_m = -0.1055613458f * a - 0.0638541728f * b;
211 float k_s = -0.0894841775f * a - 1.2914855480f * b;
212
213 float l_dt = dL + dC * k_l;
214 float m_dt = dL + dC * k_m;
215 float s_dt = dL + dC * k_s;
216
217
218 // If higher accuracy is required, 2 or 3 iterations of the following block can be used:
219 {
220 float L = L0 * (1.f - t) + t * L1;
221 float C = t * C1;
222
223 float l_ = L + C * k_l;
224 float m_ = L + C * k_m;
225 float s_ = L + C * k_s;
226
227 float l = l_ * l_ * l_;
228 float m = m_ * m_ * m_;
229 float s = s_ * s_ * s_;
230
231 float ldt = 3 * l_dt * l_ * l_;
232 float mdt = 3 * m_dt * m_ * m_;
233 float sdt = 3 * s_dt * s_ * s_;
234
235 float ldt2 = 6 * l_dt * l_dt * l_;
236 float mdt2 = 6 * m_dt * m_dt * m_;
237 float sdt2 = 6 * s_dt * s_dt * s_;
238
239 float r = 4.0767416621f * l - 3.3077115913f * m + 0.2309699292f * s - 1;
240 float r1 = 4.0767416621f * ldt - 3.3077115913f * mdt + 0.2309699292f * sdt;
241 float r2 = 4.0767416621f * ldt2 - 3.3077115913f * mdt2 + 0.2309699292f * sdt2;
242
243 float u_r = r1 / (r1 * r1 - 0.5f * r * r2);
244 float t_r = -r * u_r;
245
246 float g = -1.2684380046f * l + 2.6097574011f * m - 0.3413193965f * s - 1;
247 float g1 = -1.2684380046f * ldt + 2.6097574011f * mdt - 0.3413193965f * sdt;
248 float g2 = -1.2684380046f * ldt2 + 2.6097574011f * mdt2 - 0.3413193965f * sdt2;
249
250 float u_g = g1 / (g1 * g1 - 0.5f * g * g2);
251 float t_g = -g * u_g;
252
253 b = -0.0041960863f * l - 0.7034186147f * m + 1.7076147010f * s - 1;
254 float b1 = -0.0041960863f * ldt - 0.7034186147f * mdt + 1.7076147010f * sdt;
255 float b2 = -0.0041960863f * ldt2 - 0.7034186147f * mdt2 + 1.7076147010f * sdt2;
256
257 float u_b = b1 / (b1 * b1 - 0.5f * b * b2);
258 float t_b = -b * u_b;
259
260 t_r = u_r >= 0.f ? t_r : FLT_MAX;
261 t_g = u_g >= 0.f ? t_g : FLT_MAX;
262 t_b = u_b >= 0.f ? t_b : FLT_MAX;
263
264 t += fmin(t_r, fmin(t_g, t_b));
265 }
266 }
267 }
268
269 return t;
270}
271
272static float find_gamut_intersection(float a, float b, float L1, float C1, float L0)
273{
274 // Find the cusp of the gamut triangle
275 LC cusp = find_cusp(a, b);
276
277 return find_gamut_intersection(a, b, L1, C1, L0, cusp);
278}
279
280static RGB gamut_clip_preserve_chroma(RGB rgb)
281{
282 if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && rgb.r > 0 && rgb.g > 0 && rgb.b > 0)
283 return rgb;
284
285 Lab lab = linear_srgb_to_oklab(rgb);
286
287 float L = lab.L;
288 float eps = 0.00001f;
289 float C = fmax(eps, sqrtf(lab.a * lab.a + lab.b * lab.b));
290 float a_ = lab.a / C;
291 float b_ = lab.b / C;
292
293 float L0 = clamp(L, 0, 1);
294
295 float t = find_gamut_intersection(a_, b_, L, C, L0);
296 float L_clipped = L0 * (1 - t) + t * L;
297 float C_clipped = t * C;
298
299 return oklab_to_linear_srgb({ L_clipped, C_clipped * a_, C_clipped * b_ });
300}
301
302static RGB gamut_clip_project_to_0_5(RGB rgb)
303{
304 if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && rgb.r > 0 && rgb.g > 0 && rgb.b > 0)
305 return rgb;
306
307 Lab lab = linear_srgb_to_oklab(rgb);
308
309 float L = lab.L;
310 float eps = 0.00001f;
311 float C = fmax(eps, sqrtf(lab.a * lab.a + lab.b * lab.b));
312 float a_ = lab.a / C;
313 float b_ = lab.b / C;
314
315 float L0 = 0.5;
316
317 float t = find_gamut_intersection(a_, b_, L, C, L0);
318 float L_clipped = L0 * (1 - t) + t * L;
319 float C_clipped = t * C;
320
321 return oklab_to_linear_srgb({ L_clipped, C_clipped * a_, C_clipped * b_ });
322}
323
324static RGB gamut_clip_project_to_L_cusp(RGB rgb)
325{
326 if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && rgb.r > 0 && rgb.g > 0 && rgb.b > 0)
327 return rgb;
328
329 Lab lab = linear_srgb_to_oklab(rgb);
330
331 float L = lab.L;
332 float eps = 0.00001f;
333 float C = fmax(eps, sqrtf(lab.a * lab.a + lab.b * lab.b));
334 float a_ = lab.a / C;
335 float b_ = lab.b / C;
336
337 // The cusp is computed here and in find_gamut_intersection, an optimized solution would only compute it once.
338 LC cusp = find_cusp(a_, b_);
339
340 float L0 = cusp.L;
341
342 float t = find_gamut_intersection(a_, b_, L, C, L0);
343
344 float L_clipped = L0 * (1 - t) + t * L;
345 float C_clipped = t * C;
346
347 return oklab_to_linear_srgb({ L_clipped, C_clipped * a_, C_clipped * b_ });
348}
349
350static RGB gamut_clip_adaptive_L0_0_5(RGB rgb, float alpha = 0.05f)
351{
352 if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && rgb.r > 0 && rgb.g > 0 && rgb.b > 0)
353 return rgb;
354
355 Lab lab = linear_srgb_to_oklab(rgb);
356
357 float L = lab.L;
358 float eps = 0.00001f;
359 float C = fmax(eps, sqrtf(lab.a * lab.a + lab.b * lab.b));
360 float a_ = lab.a / C;
361 float b_ = lab.b / C;
362
363 float Ld = L - 0.5f;
364 float e1 = 0.5f + fabs(Ld) + alpha * C;
365 float L0 = 0.5f * (1.f + sgn(Ld) * (e1 - sqrtf(e1 * e1 - 2.f * fabs(Ld))));
366
367 float t = find_gamut_intersection(a_, b_, L, C, L0);
368 float L_clipped = L0 * (1.f - t) + t * L;
369 float C_clipped = t * C;
370
371 return oklab_to_linear_srgb({ L_clipped, C_clipped * a_, C_clipped * b_ });
372}
373
374static RGB gamut_clip_adaptive_L0_L_cusp(RGB rgb, float alpha = 0.05f)
375{
376 if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && rgb.r > 0 && rgb.g > 0 && rgb.b > 0)
377 return rgb;
378
379 Lab lab = linear_srgb_to_oklab(rgb);
380
381 float L = lab.L;
382 float eps = 0.00001f;
383 float C = fmax(eps, sqrtf(lab.a * lab.a + lab.b * lab.b));
384 float a_ = lab.a / C;
385 float b_ = lab.b / C;
386
387 // The cusp is computed here and in find_gamut_intersection, an optimized solution would only compute it once.
388 LC cusp = find_cusp(a_, b_);
389
390 float Ld = L - cusp.L;
391 float k = 2.f * (Ld > 0 ? 1.f - cusp.L : cusp.L);
392
393 float e1 = 0.5f * k + fabs(Ld) + alpha * C / k;
394 float L0 = cusp.L + 0.5f * (sgn(Ld) * (e1 - sqrtf(e1 * e1 - 2.f * k * fabs(Ld))));
395
396 float t = find_gamut_intersection(a_, b_, L, C, L0);
397 float L_clipped = L0 * (1.f - t) + t * L;
398 float C_clipped = t * C;
399
400 return oklab_to_linear_srgb({ L_clipped, C_clipped * a_, C_clipped * b_ });
401}
402
403static float toe(float x)
404{
405 constexpr float k_1 = 0.206f;
406 constexpr float k_2 = 0.03f;
407 constexpr float k_3 = (1.f + k_1) / (1.f + k_2);
408 return 0.5f * (k_3 * x - k_1 + sqrtf((k_3 * x - k_1) * (k_3 * x - k_1) + 4 * k_2 * k_3 * x));
409}
410
411static float toe_inv(float x)
412{
413 constexpr float k_1 = 0.206f;
414 constexpr float k_2 = 0.03f;
415 constexpr float k_3 = (1.f + k_1) / (1.f + k_2);
416 return (x * x + k_1 * x) / (k_3 * (x + k_2));
417}
418
419static ST to_ST(LC cusp)
420{
421 float L = cusp.L;
422 float C = cusp.C;
423 return { C / L, C / (1 - L) };
424}
425
426// Returns a smooth approximation of the location of the cusp
427// This polynomial was created by an optimization process
428// It has been designed so that S_mid < S_max and T_mid < T_max
429static ST get_ST_mid(float a_, float b_)
430{
431 float S = 0.11516993f + 1.f / (
432 +7.44778970f + 4.15901240f * b_
433 + a_ * (-2.19557347f + 1.75198401f * b_
434 + a_ * (-2.13704948f - 10.02301043f * b_
435 + a_ * (-4.24894561f + 5.38770819f * b_ + 4.69891013f * a_
436 )))
437 );
438
439 float T = 0.11239642f + 1.f / (
440 +1.61320320f - 0.68124379f * b_
441 + a_ * (+0.40370612f + 0.90148123f * b_
442 + a_ * (-0.27087943f + 0.61223990f * b_
443 + a_ * (+0.00299215f - 0.45399568f * b_ - 0.14661872f * a_
444 )))
445 );
446
447 return { S, T };
448}
449
450struct Cs { float C_0; float C_mid; float C_max; };
451static Cs get_Cs(float L, float a_, float b_)
452{
453 LC cusp = find_cusp(a_, b_);
454
455 float C_max = find_gamut_intersection(a_, b_, L, 1, L, cusp);
456 ST ST_max = to_ST(cusp);
457
458 // Scale factor to compensate for the curved part of gamut shape:
459 float k = C_max / fmin((L * ST_max.S), (1 - L) * ST_max.T);
460
461 float C_mid;
462 {
463 ST ST_mid = get_ST_mid(a_, b_);
464
465 // Use a soft minimum function, instead of a sharp triangle shape to get a smooth value for chroma.
466 float C_a = L * ST_mid.S;
467 float C_b = (1.f - L) * ST_mid.T;
468 C_mid = 0.9f * k * sqrtf(sqrtf(1.f / (1.f / (C_a * C_a * C_a * C_a) + 1.f / (C_b * C_b * C_b * C_b))));
469 }
470
471 float C_0;
472 {
473 // for C_0, the shape is independent of hue, so ST are constant. Values picked to roughly be the average values of ST.
474 float C_a = L * 0.4f;
475 float C_b = (1.f - L) * 0.8f;
476
477 // Use a soft minimum function, instead of a sharp triangle shape to get a smooth value for chroma.
478 C_0 = sqrtf(1.f / (1.f / (C_a * C_a) + 1.f / (C_b * C_b)));
479 }
480
481 return { C_0, C_mid, C_max };
482}
483
484static RGB okhsl_to_srgb(HSL hsl)
485{
486 float h = hsl.h;
487 float s = hsl.s;
488 float l = hsl.l;
489
490 if (l == 1.0f)
491 {
492 return { 1.f, 1.f, 1.f };
493 }
494
495 else if (l == 0.f)
496 {
497 return { 0.f, 0.f, 0.f };
498 }
499
500 float a_ = cosf(2.f * pi * h);
501 float b_ = sinf(2.f * pi * h);
502 float L = toe_inv(l);
503
504 Cs cs = get_Cs(L, a_, b_);
505 float C_0 = cs.C_0;
506 float C_mid = cs.C_mid;
507 float C_max = cs.C_max;
508
509 float mid = 0.8f;
510 float mid_inv = 1.25f;
511
512 float C, t, k_0, k_1, k_2;
513
514 if (s < mid)
515 {
516 t = mid_inv * s;
517
518 k_1 = mid * C_0;
519 k_2 = (1.f - k_1 / C_mid);
520
521 C = t * k_1 / (1.f - k_2 * t);
522 }
523 else
524 {
525 t = (s - mid)/ (1 - mid);
526
527 k_0 = C_mid;
528 k_1 = (1.f - mid) * C_mid * C_mid * mid_inv * mid_inv / C_0;
529 k_2 = (1.f - (k_1) / (C_max - C_mid));
530
531 C = k_0 + t * k_1 / (1.f - k_2 * t);
532 }
533
534 RGB rgb = oklab_to_linear_srgb({ L, C * a_, C * b_ });
535 return {
536 srgb_transfer_function(rgb.r),
537 srgb_transfer_function(rgb.g),
538 srgb_transfer_function(rgb.b),
539 };
540}
541
542static HSL srgb_to_okhsl(RGB rgb)
543{
544 Lab lab = linear_srgb_to_oklab({
545 srgb_transfer_function_inv(rgb.r),
546 srgb_transfer_function_inv(rgb.g),
547 srgb_transfer_function_inv(rgb.b)
548 });
549
550 float C = sqrtf(lab.a * lab.a + lab.b * lab.b);
551 float a_ = lab.a / C;
552 float b_ = lab.b / C;
553
554 float L = lab.L;
555 float h = 0.5f + 0.5f * atan2f(-lab.b, -lab.a) / pi;
556
557 Cs cs = get_Cs(L, a_, b_);
558 float C_0 = cs.C_0;
559 float C_mid = cs.C_mid;
560 float C_max = cs.C_max;
561
562 // Inverse of the interpolation in okhsl_to_srgb:
563
564 float mid = 0.8f;
565 float mid_inv = 1.25f;
566
567 float s;
568 if (C < C_mid)
569 {
570 float k_1 = mid * C_0;
571 float k_2 = (1.f - k_1 / C_mid);
572
573 float t = C / (k_1 + k_2 * C);
574 s = t * mid;
575 }
576 else
577 {
578 float k_0 = C_mid;
579 float k_1 = (1.f - mid) * C_mid * C_mid * mid_inv * mid_inv / C_0;
580 float k_2 = (1.f - (k_1) / (C_max - C_mid));
581
582 float t = (C - k_0) / (k_1 + k_2 * (C - k_0));
583 s = mid + (1.f - mid) * t;
584 }
585
586 float l = toe(L);
587 return { h, s, l };
588}
589
590
591static RGB okhsv_to_srgb(HSV hsv)
592{
593 float h = hsv.h;
594 float s = hsv.s;
595 float v = hsv.v;
596
597 float a_ = cosf(2.f * pi * h);
598 float b_ = sinf(2.f * pi * h);
599
600 LC cusp = find_cusp(a_, b_);
601 ST ST_max = to_ST(cusp);
602 float S_max = ST_max.S;
603 float T_max = ST_max.T;
604 float S_0 = 0.5f;
605 float k = 1 - S_0 / S_max;
606
607 // first we compute L and V as if the gamut is a perfect triangle:
608
609 // L, C when v==1:
610 float L_v = 1 - s * S_0 / (S_0 + T_max - T_max * k * s);
611 float C_v = s * T_max * S_0 / (S_0 + T_max - T_max * k * s);
612
613 float L = v * L_v;
614 float C = v * C_v;
615
616 // then we compensate for both toe and the curved top part of the triangle:
617 float L_vt = toe_inv(L_v);
618 float C_vt = C_v * L_vt / L_v;
619
620 float L_new = toe_inv(L);
621 C = C * L_new / L;
622 L = L_new;
623
624 RGB rgb_scale = oklab_to_linear_srgb({ L_vt, a_ * C_vt, b_ * C_vt });
625 float scale_L = cbrtf(1.f / fmax(fmax(rgb_scale.r, rgb_scale.g), fmax(rgb_scale.b, 0.f)));
626
627 L = L * scale_L;
628 C = C * scale_L;
629
630 RGB rgb = oklab_to_linear_srgb({ L, C * a_, C * b_ });
631 return {
632 srgb_transfer_function(rgb.r),
633 srgb_transfer_function(rgb.g),
634 srgb_transfer_function(rgb.b),
635 };
636}
637
638static HSV srgb_to_okhsv(RGB rgb)
639{
640 Lab lab = linear_srgb_to_oklab({
641 srgb_transfer_function_inv(rgb.r),
642 srgb_transfer_function_inv(rgb.g),
643 srgb_transfer_function_inv(rgb.b)
644 });
645
646 float C = sqrtf(lab.a * lab.a + lab.b * lab.b);
647 float a_ = lab.a / C;
648 float b_ = lab.b / C;
649
650 float L = lab.L;
651 float h = 0.5f + 0.5f * atan2f(-lab.b, -lab.a) / pi;
652
653 LC cusp = find_cusp(a_, b_);
654 ST ST_max = to_ST(cusp);
655 float S_max = ST_max.S;
656 float T_max = ST_max.T;
657 float S_0 = 0.5f;
658 float k = 1 - S_0 / S_max;
659
660 // first we find L_v, C_v, L_vt and C_vt
661
662 float t = T_max / (C + L * T_max);
663 float L_v = t * L;
664 float C_v = t * C;
665
666 float L_vt = toe_inv(L_v);
667 float C_vt = C_v * L_vt / L_v;
668
669 // we can then use these to invert the step that compensates for the toe and the curved top part of the triangle:
670 RGB rgb_scale = oklab_to_linear_srgb({ L_vt, a_ * C_vt, b_ * C_vt });
671 float scale_L = cbrtf(1.f / fmax(fmax(rgb_scale.r, rgb_scale.g), fmax(rgb_scale.b, 0.f)));
672
673 L = L / scale_L;
674 C = C / scale_L;
675
676 C = C * toe(L) / L;
677 L = toe(L);
678
679 // we can now compute v and s:
680
681 float v = L / L_v;
682 float s = (S_0 + T_max) * C_v / ((T_max * S_0) + T_max * k * C_v);
683
684 return { h, s, v };
685}
686
687};
688#endif // OK_COLOR_H
689