LSST Applications 24.1.6,g063fba187b+56b85ce14a,g0f08755f38+df8a265115,g12f32b3c4e+891a09f10d,g1524ad2192+7a5d7b3fbd,g1653933729+a8ce1bb630,g168dd56ebc+a8ce1bb630,g28da252d5a+07cb1400be,g2bbee38e9b+ae03bbfc84,g2bc492864f+ae03bbfc84,g3156d2b45e+6e55a43351,g347aa1857d+ae03bbfc84,g35bb328faa+a8ce1bb630,g3a166c0a6a+ae03bbfc84,g3e281a1b8c+c5dd892a6c,g414038480c+6b9177ef31,g41af890bb2+8f257c4c0b,g781aacb6e4+a8ce1bb630,g7af13505b9+7137b3b17d,g80478fca09+6df6903293,g82479be7b0+091ce1d07f,g858d7b2824+df8a265115,g89c8672015+f4add4ffd5,g9125e01d80+a8ce1bb630,g9726552aa6+414189b318,ga5288a1d22+4a2bca08d7,gacef1a1666+c9a8ff65f4,gb58c049af0+d64f4d3760,gbcfae0f0a0+de1d42d831,gc28159a63d+ae03bbfc84,gcf0d15dbbd+72117bf34e,gda6a2b7d83+72117bf34e,gdaeeff99f8+1711a396fd,ge500cccec5+c8c9c9af63,ge79ae78c31+ae03bbfc84,gf0baf85859+c1f95f4921,gfa517265be+df8a265115,gfa999e8aa5+17cd334064,gfb92a5be7c+df8a265115
LSST Data Management Base Package
Loading...
Searching...
No Matches
Classes | Functions
lsst::cpputils::details Namespace Reference

Classes

struct  Lab
 
struct  LC
 
struct  RGB
 

Functions

Lab linear_srgb_to_oklab (RGB c)
 
Lab linear_displayP3_to_oklab (RGB c)
 
RGB oklab_to_linear_srgb (Lab c)
 
RGB oklab_to_linear_displayP3 (Lab c)
 
float compute_max_saturation (float a, float b)
 
LC find_cusp (float a, float b)
 
float find_gamut_intersection (float a, float b, float L1, float C1, float L0)
 
float clamp (float x, float min, float max)
 
float sgn (float x)
 
RGB gamut_clip_preserve_chroma (RGB rgb)
 
RGB gamut_clip_project_to_0_5 (RGB rgb)
 
RGB gamut_clip_project_to_L_cusp (RGB rgb)
 
RGB gamut_clip_adaptive_L0_0_5 (RGB rgb, float alpha=0.05f)
 
RGB gamut_clip_adaptive_L0_L_cusp (RGB rgb, float alpha=0.05f)
 

Function Documentation

◆ clamp()

float lsst::cpputils::details::clamp ( float x,
float min,
float max )

Definition at line 281 of file _oklabTools.h.

282{
283 if (x < min)
284 return min;
285 if (x > max)
286 return max;
287
288 return x;
289}
int min
int max

◆ compute_max_saturation()

float lsst::cpputils::details::compute_max_saturation ( float a,
float b )

Definition at line 109 of file _oklabTools.h.

110{
111 // Max saturation will be when one of r, g or b goes below zero.
112
113 // Select different coefficients depending on which component goes below zero first
114 float k0, k1, k2, k3, k4, wl, wm, ws;
115
116 if (-1.88170328f * a - 0.80936493f * b > 1)
117 {
118 // Red component
119 k0 = +1.19086277f; k1 = +1.76576728f; k2 = +0.59662641f; k3 = +0.75515197f; k4 = +0.56771245f;
120 wl = +4.0767416621f; wm = -3.3077115913f; ws = +0.2309699292f;
121 }
122 else if (1.81444104f * a - 1.19445276f * b > 1)
123 {
124 // Green component
125 k0 = +0.73956515f; k1 = -0.45954404f; k2 = +0.08285427f; k3 = +0.12541070f; k4 = +0.14503204f;
126 wl = -1.2684380046f; wm = +2.6097574011f; ws = -0.3413193965f;
127 }
128 else
129 {
130 // Blue component
131 k0 = +1.35733652f; k1 = -0.00915799f; k2 = -1.15130210f; k3 = -0.50559606f; k4 = +0.00692167f;
132 wl = -0.0041960863f; wm = -0.7034186147f; ws = +1.7076147010f;
133 }
134
135 // Approximate max saturation using a polynomial:
136 float S = k0 + k1 * a + k2 * b + k3 * a * a + k4 * a * b;
137
138 // Do one step Halley's method to get closer
139 // this gives an error less than 10e6, except for some blue hues where the dS/dh is close to infinite
140 // this should be sufficient for most applications, otherwise do two/three steps
141
142 float k_l = +0.3963377774f * a + 0.2158037573f * b;
143 float k_m = -0.1055613458f * a - 0.0638541728f * b;
144 float k_s = -0.0894841775f * a - 1.2914855480f * b;
145
146 {
147 float l_ = 1.f + S * k_l;
148 float m_ = 1.f + S * k_m;
149 float s_ = 1.f + S * k_s;
150
151 float l = l_ * l_ * l_;
152 float m = m_ * m_ * m_;
153 float s = s_ * s_ * s_;
154
155 float l_dS = 3.f * k_l * l_ * l_;
156 float m_dS = 3.f * k_m * m_ * m_;
157 float s_dS = 3.f * k_s * s_ * s_;
158
159 float l_dS2 = 6.f * k_l * k_l * l_;
160 float m_dS2 = 6.f * k_m * k_m * m_;
161 float s_dS2 = 6.f * k_s * k_s * s_;
162
163 float f = wl * l + wm * m + ws * s;
164 float f1 = wl * l_dS + wm * m_dS + ws * s_dS;
165 float f2 = wl * l_dS2 + wm * m_dS2 + ws * s_dS2;
166
167 S = S - f * f1 / (f1*f1 - 0.5f * f * f2);
168 }
169
170 return S;
171}
int m
Definition SpanSet.cc:48
table::Key< int > b

◆ find_cusp()

LC lsst::cpputils::details::find_cusp ( float a,
float b )

Definition at line 176 of file _oklabTools.h.

177{
178 // First, find the maximum saturation (saturation S = C/L)
179 float S_cusp = compute_max_saturation(a, b);
180
181 // Convert to linear sRGB to find the first point where at least one of r,g or b >= 1:
182 RGB rgb_at_max = oklab_to_linear_displayP3({ 1, S_cusp * a, S_cusp * b });
183 float L_cusp = cbrtf(1.f / fmax(fmax(rgb_at_max.r, rgb_at_max.g), rgb_at_max.b));
184 float C_cusp = L_cusp * S_cusp;
185
186 return { L_cusp , C_cusp };
187}
table::Key< int > a
RGB oklab_to_linear_displayP3(Lab c)
Definition _oklabTools.h:89
float compute_max_saturation(float a, float b)

◆ find_gamut_intersection()

float lsst::cpputils::details::find_gamut_intersection ( float a,
float b,
float L1,
float C1,
float L0 )

Definition at line 193 of file _oklabTools.h.

194{
195 // Find the cusp of the gamut triangle
196 LC cusp = find_cusp(a, b);
197
198 // Find the intersection for upper and lower half seprately
199 float t;
200 if (((L1 - L0) * cusp.C - (cusp.L - L0) * C1) <= 0.f)
201 {
202 // Lower half
203
204 t = cusp.C * L0 / (C1 * cusp.L + cusp.C * (L0 - L1));
205 }
206 else
207 {
208 // Upper half
209
210 // First intersect with triangle
211 t = cusp.C * (L0 - 1.f) / (C1 * (cusp.L - 1.f) + cusp.C * (L0 - L1));
212
213 // Then one step Halley's method
214 {
215 float dL = L1 - L0;
216 float dC = C1;
217
218 float k_l = +0.3963377774f * a + 0.2158037573f * b;
219 float k_m = -0.1055613458f * a - 0.0638541728f * b;
220 float k_s = -0.0894841775f * a - 1.2914855480f * b;
221
222 float l_dt = dL + dC * k_l;
223 float m_dt = dL + dC * k_m;
224 float s_dt = dL + dC * k_s;
225
226
227 // If higher accuracy is required, 2 or 3 iterations of the following block can be used:
228 {
229 float L = L0 * (1.f - t) + t * L1;
230 float C = t * C1;
231
232 float l_ = L + C * k_l;
233 float m_ = L + C * k_m;
234 float s_ = L + C * k_s;
235
236 float l = l_ * l_ * l_;
237 float m = m_ * m_ * m_;
238 float s = s_ * s_ * s_;
239
240 float ldt = 3 * l_dt * l_ * l_;
241 float mdt = 3 * m_dt * m_ * m_;
242 float sdt = 3 * s_dt * s_ * s_;
243
244 float ldt2 = 6 * l_dt * l_dt * l_;
245 float mdt2 = 6 * m_dt * m_dt * m_;
246 float sdt2 = 6 * s_dt * s_dt * s_;
247
248 float r = 4.0767416621f * l - 3.3077115913f * m + 0.2309699292f * s - 1;
249 float r1 = 4.0767416621f * ldt - 3.3077115913f * mdt + 0.2309699292f * sdt;
250 float r2 = 4.0767416621f * ldt2 - 3.3077115913f * mdt2 + 0.2309699292f * sdt2;
251
252 float u_r = r1 / (r1 * r1 - 0.5f * r * r2);
253 float t_r = -r * u_r;
254
255 float g = -1.2684380046f * l + 2.6097574011f * m - 0.3413193965f * s - 1;
256 float g1 = -1.2684380046f * ldt + 2.6097574011f * mdt - 0.3413193965f * sdt;
257 float g2 = -1.2684380046f * ldt2 + 2.6097574011f * mdt2 - 0.3413193965f * sdt2;
258
259 float u_g = g1 / (g1 * g1 - 0.5f * g * g2);
260 float t_g = -g * u_g;
261
262 float b = -0.0041960863f * l - 0.7034186147f * m + 1.7076147010f * s - 1;
263 float b1 = -0.0041960863f * ldt - 0.7034186147f * mdt + 1.7076147010f * sdt;
264 float b2 = -0.0041960863f * ldt2 - 0.7034186147f * mdt2 + 1.7076147010f * sdt2;
265
266 float u_b = b1 / (b1 * b1 - 0.5f * b * b2);
267 float t_b = -b * u_b;
268
269 t_r = u_r >= 0.f ? t_r : FLT_MAX;
270 t_g = u_g >= 0.f ? t_g : FLT_MAX;
271 t_b = u_b >= 0.f ? t_b : FLT_MAX;
272
273 t += fmin(t_r, fmin(t_g, t_b));
274 }
275 }
276 }
277
278 return t;
279}
LC find_cusp(float a, float b)

◆ gamut_clip_adaptive_L0_0_5()

RGB lsst::cpputils::details::gamut_clip_adaptive_L0_0_5 ( RGB rgb,
float alpha = 0.05f )

Definition at line 366 of file _oklabTools.h.

367{
368 if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && rgb.r > 0 && rgb.g > 0 && rgb.b > 0)
369 return rgb;
370
372
373 float L = lab.L;
374 float eps = 0.00001f;
375 float C = fmax(eps, sqrtf(lab.a * lab.a + lab.b * lab.b));
376 float a_ = lab.a / C;
377 float b_ = lab.b / C;
378
379 float Ld = L - 0.5f;
380 float e1 = 0.5f + fabs(Ld) + alpha * C;
381 float L0 = 0.5f*(1.f + sgn(Ld)*(e1 - sqrtf(e1*e1 - 2.f *fabs(Ld))));
382
383 float t = find_gamut_intersection(a_, b_, L, C, L0);
384 float L_clipped = L0 * (1.f - t) + t * L;
385 float C_clipped = t * C;
386
387 return oklab_to_linear_displayP3({ L_clipped, C_clipped * a_, C_clipped * b_ });
388}
Lab linear_displayP3_to_oklab(RGB c)
Definition _oklabTools.h:54
float find_gamut_intersection(float a, float b, float L1, float C1, float L0)

◆ gamut_clip_adaptive_L0_L_cusp()

RGB lsst::cpputils::details::gamut_clip_adaptive_L0_L_cusp ( RGB rgb,
float alpha = 0.05f )

Definition at line 390 of file _oklabTools.h.

391{
392 if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && rgb.r > 0 && rgb.g > 0 && rgb.b > 0)
393 return rgb;
394
396
397 float L = lab.L;
398 float eps = 0.00001f;
399 float C = fmax(eps, sqrtf(lab.a * lab.a + lab.b * lab.b));
400 float a_ = lab.a / C;
401 float b_ = lab.b / C;
402
403 // The cusp is computed here and in find_gamut_intersection, an optimized solution would only compute it once.
404 LC cusp = find_cusp(a_, b_);
405
406 float Ld = L - cusp.L;
407 float k = 2.f * (Ld > 0 ? 1.f - cusp.L : cusp.L);
408
409 float e1 = 0.5f*k + fabs(Ld) + alpha * C/k;
410 float L0 = cusp.L + 0.5f * (sgn(Ld) * (e1 - sqrtf(e1 * e1 - 2.f * k * fabs(Ld))));
411
412 float t = find_gamut_intersection(a_, b_, L, C, L0);
413 float L_clipped = L0 * (1.f - t) + t * L;
414 float C_clipped = t * C;
415
416 return oklab_to_linear_displayP3({ L_clipped, C_clipped * a_, C_clipped * b_ });
417}

◆ gamut_clip_preserve_chroma()

RGB lsst::cpputils::details::gamut_clip_preserve_chroma ( RGB rgb)

Definition at line 296 of file _oklabTools.h.

297{
298 if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && rgb.r > 0 && rgb.g > 0 && rgb.b > 0)
299 return rgb;
300
302
303 float L = lab.L;
304 float eps = 0.00001f;
305 float C = fmax(eps, sqrtf(lab.a * lab.a + lab.b * lab.b));
306 float a_ = lab.a / C;
307 float b_ = lab.b / C;
308
309 float L0 = clamp(L, 0, 1);
310
311 float t = find_gamut_intersection(a_, b_, L, C, L0);
312 float L_clipped = L0 * (1 - t) + t * L;
313 float C_clipped = t * C;
314
315 return oklab_to_linear_displayP3({ L_clipped, C_clipped * a_, C_clipped * b_ });
316}
float clamp(float x, float min, float max)

◆ gamut_clip_project_to_0_5()

RGB lsst::cpputils::details::gamut_clip_project_to_0_5 ( RGB rgb)

Definition at line 318 of file _oklabTools.h.

319{
320 if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && rgb.r > 0 && rgb.g > 0 && rgb.b > 0)
321 return rgb;
322
324
325 float L = lab.L;
326 float eps = 0.00001f;
327 float C = fmax(eps, sqrtf(lab.a * lab.a + lab.b * lab.b));
328 float a_ = lab.a / C;
329 float b_ = lab.b / C;
330
331 float L0 = 0.5;
332
333 float t = find_gamut_intersection(a_, b_, L, C, L0);
334 float L_clipped = L0 * (1 - t) + t * L;
335 float C_clipped = t * C;
336
337 return oklab_to_linear_displayP3({ L_clipped, C_clipped * a_, C_clipped * b_ });
338}

◆ gamut_clip_project_to_L_cusp()

RGB lsst::cpputils::details::gamut_clip_project_to_L_cusp ( RGB rgb)

Definition at line 340 of file _oklabTools.h.

341{
342 if (rgb.r < 1 && rgb.g < 1 && rgb.b < 1 && rgb.r > 0 && rgb.g > 0 && rgb.b > 0)
343 return rgb;
344
346
347 float L = lab.L;
348 float eps = 0.00001f;
349 float C = fmax(eps, sqrtf(lab.a * lab.a + lab.b * lab.b));
350 float a_ = lab.a / C;
351 float b_ = lab.b / C;
352
353 // The cusp is computed here and in find_gamut_intersection, an optimized solution would only compute it once.
354 LC cusp = find_cusp(a_, b_);
355
356 float L0 = cusp.L;
357
358 float t = find_gamut_intersection(a_, b_, L, C, L0);
359
360 float L_clipped = L0 * (1 - t) + t * L;
361 float C_clipped = t * C;
362
363 return oklab_to_linear_displayP3({ L_clipped, C_clipped * a_, C_clipped * b_ });
364}

◆ linear_displayP3_to_oklab()

Lab lsst::cpputils::details::linear_displayP3_to_oklab ( RGB c)

Definition at line 54 of file _oklabTools.h.

55{
56 float l = 0.48132729f * c.r + 0.46206791f * c.g + 00.0564956f * c.b;
57 float m = 0.2288381f * c.r + 0.6532344f * c.g + 0.11795441f * c.b;
58 float s = 0.08398602f * c.r + 0.22427279f * c.g + 0.69222084f * c.b;
59
60 float l_ = cbrtf(l);
61 float m_ = cbrtf(m);
62 float s_ = cbrtf(s);
63
64 return {
65 0.2104542553f * l_ + 0.7936177850f * m_ - 0.0040720468f * s_,
66 1.9779984951f * l_ - 2.4285922050f * m_ + 0.4505937099f * s_,
67 0.0259040371f * l_ + 0.7827717662f * m_ - 0.8086757660f * s_,
68 };
69}

◆ linear_srgb_to_oklab()

Lab lsst::cpputils::details::linear_srgb_to_oklab ( RGB c)

Definition at line 37 of file _oklabTools.h.

38{
39 float l = 0.4122214708f * c.r + 0.5363325363f * c.g + 0.0514459929f * c.b;
40 float m = 0.2119034982f * c.r + 0.6806995451f * c.g + 0.1073969566f * c.b;
41 float s = 0.0883024619f * c.r + 0.2817188376f * c.g + 0.6299787005f * c.b;
42
43 float l_ = cbrtf(l);
44 float m_ = cbrtf(m);
45 float s_ = cbrtf(s);
46
47 return {
48 0.2104542553f * l_ + 0.7936177850f * m_ - 0.0040720468f * s_,
49 1.9779984951f * l_ - 2.4285922050f * m_ + 0.4505937099f * s_,
50 0.0259040371f * l_ + 0.7827717662f * m_ - 0.8086757660f * s_,
51 };
52}

◆ oklab_to_linear_displayP3()

RGB lsst::cpputils::details::oklab_to_linear_displayP3 ( Lab c)

Definition at line 89 of file _oklabTools.h.

90{
91 float l_ = c.L + 0.3963377774f * c.a + 0.2158037573f * c.b;
92 float m_ = c.L - 0.1055613458f * c.a - 0.0638541728f * c.b;
93 float s_ = c.L - 0.0894841775f * c.a - 1.2914855480f * c.b;
94
95 float l = l_ * l_ * l_;
96 float m = m_ * m_ * m_;
97 float s = s_ * s_ * s_;
98
99 return {
100 3.12811053f * l - 2.25707502f * m + 0.12930479f * s,
101 -1.09112816f * l + 2.41326676f * m - 0.32216817f * s,
102 -0.02601365f * l - 0.50802765f * m + 1.53331668f * s,
103 };
104}

◆ oklab_to_linear_srgb()

RGB lsst::cpputils::details::oklab_to_linear_srgb ( Lab c)

Definition at line 72 of file _oklabTools.h.

73{
74 float l_ = c.L + 0.3963377774f * c.a + 0.2158037573f * c.b;
75 float m_ = c.L - 0.1055613458f * c.a - 0.0638541728f * c.b;
76 float s_ = c.L - 0.0894841775f * c.a - 1.2914855480f * c.b;
77
78 float l = l_ * l_ * l_;
79 float m = m_ * m_ * m_;
80 float s = s_ * s_ * s_;
81
82 return {
83 +4.0767416621f * l - 3.3077115913f * m + 0.2309699292f * s,
84 -1.2684380046f * l + 2.6097574011f * m - 0.3413193965f * s,
85 -0.0041960863f * l - 0.7034186147f * m + 1.7076147010f * s,
86 };
87}

◆ sgn()

float lsst::cpputils::details::sgn ( float x)

Definition at line 291 of file _oklabTools.h.

292{
293 return (float)(0.f < x) - (float)(x < 0.f);
294}