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