?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Bug 1799258 - Support outByIn.size()<2 in SampleOutByIn. r=bradwerth
[gecko.git] / gfx / gl / Colorspaces.h
blob8f36854d2d4f3ab277b597401a723638fa5b676d
1 /* This Source Code Form is subject to the terms of the Mozilla Public
2 * License, v. 2.0. If a copy of the MPL was not distributed with this
3 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
4
5 #ifndef MOZILLA_GFX_GL_COLORSPACES_H_
6 #define MOZILLA_GFX_GL_COLORSPACES_H_
7
8 // Reference: https://hackmd.io/0wkiLmP7RWOFjcD13M870A
9
10 // We are going to be doing so, so many transforms, so descriptive labels are
11 // critical.
12
13 // Colorspace background info: https://hackmd.io/0wkiLmP7RWOFjcD13M870A
14
15 #include <algorithm>
16 #include <array>
17 #include <cmath>
18 #include <cstdint>
19 #include <cstdlib>
20 #include <functional>
21 #include <optional>
22 #include <vector>
23
24 #include "AutoMappable.h"
25 #include "mozilla/Assertions.h"
26 #include "mozilla/Attributes.h"
27 #include "mozilla/Span.h"
28
29 #ifdef DEBUG
30 # define ASSERT(EXPR) \
31 do { \
32 if (!(EXPR)) { \
33 __builtin_trap(); \
34 } \
35 } while (false)
36 #else
37 # define ASSERT(EXPR) (void)(EXPR)
38 #endif
39
40 struct _qcms_profile;
41 typedef struct _qcms_profile qcms_profile;
42
43 namespace mozilla::color {
44
45 struct YuvLumaCoeffs final {
46 float r = 0.2126;
47 float g = 0.7152;
48 float b = 0.0722;
49
50 auto Members() const { return std::tie(r, g, b); }
51 INLINE_AUTO_MAPPABLE(YuvLumaCoeffs)
52
53 static constexpr auto Rec709() { return YuvLumaCoeffs(); }
54
55 static constexpr auto Rec2020() {
56 return YuvLumaCoeffs{0.2627, 0.6780, 0.0593};
57 }
58 };
59
60 struct PiecewiseGammaDesc final {
61 // tf = { k * linear | linear < b
62 // { a * pow(linear, 1/g) - (1-a) | linear >= b
63
64 // Default to Srgb
65 float a = 1.055;
66 float b = 0.04045 / 12.92;
67 float g = 2.4;
68 float k = 12.92;
69
70 auto Members() const { return std::tie(a, b, g, k); }
71 INLINE_AUTO_MAPPABLE(PiecewiseGammaDesc)
72
73 static constexpr auto Srgb() { return PiecewiseGammaDesc(); }
74 static constexpr auto DisplayP3() { return Srgb(); }
75
76 static constexpr auto Rec709() {
77 return PiecewiseGammaDesc{
78 1.099,
79 0.018,
80 1.0 / 0.45, // ~2.222
81 4.5,
82 };
83 }
84 // FYI: static constexpr auto Rec2020_10bit() { return Rec709(); }
85 static constexpr auto Rec2020_12bit() {
86 return PiecewiseGammaDesc{
87 1.0993,
88 0.0181,
89 1.0 / 0.45, // ~2.222
90 4.5,
91 };
92 }
93 };
94
95 struct YcbcrDesc final {
96 float y0 = 16 / 255.0;
97 float y1 = 235 / 255.0;
98 float u0 = 128 / 255.0;
99 float uPlusHalf = 240 / 255.0;
100
101 auto Members() const { return std::tie(y0, y1, u0, uPlusHalf); }
102 INLINE_AUTO_MAPPABLE(YcbcrDesc)
103
104 static constexpr auto Narrow8() { // AKA limited/studio/tv
105 return YcbcrDesc();
106 }
107 static constexpr auto Full8() { // AKA pc
108 return YcbcrDesc{
109 0 / 255.0,
110 255 / 255.0,
111 128 / 255.0,
112 254 / 255.0,
113 };
114 }
115 static constexpr auto Float() { // Best for a LUT
116 return YcbcrDesc{0.0, 1.0, 0.5, 1.0};
117 }
118 };
119
120 struct Chromaticities final {
121 float rx = 0.640;
122 float ry = 0.330;
123 float gx = 0.300;
124 float gy = 0.600;
125 float bx = 0.150;
126 float by = 0.060;
127 // D65:
128 static constexpr float wx = 0.3127;
129 static constexpr float wy = 0.3290;
130
131 auto Members() const { return std::tie(rx, ry, gx, gy, bx, by); }
132 INLINE_AUTO_MAPPABLE(Chromaticities)
133
134 // -
135
136 static constexpr auto Rec709() { // AKA limited/studio/tv
137 return Chromaticities();
138 }
139 static constexpr auto Srgb() { return Rec709(); }
140
141 static constexpr auto Rec601_625_Pal() {
142 auto ret = Rec709();
143 ret.gx = 0.290;
144 return ret;
145 }
146 static constexpr auto Rec601_525_Ntsc() {
147 return Chromaticities{
148 0.630, 0.340, // r
149 0.310, 0.595, // g
150 0.155, 0.070, // b
151 };
152 }
153 static constexpr auto Rec2020() {
154 return Chromaticities{
155 0.708, 0.292, // r
156 0.170, 0.797, // g
157 0.131, 0.046, // b
158 };
159 }
160 static constexpr auto DisplayP3() {
161 return Chromaticities{
162 0.680, 0.320, // r
163 0.265, 0.690, // g
164 0.150, 0.060, // b
165 };
166 }
167 };
168
169 // -
170
171 struct YuvDesc final {
172 YuvLumaCoeffs yCoeffs;
173 YcbcrDesc ycbcr;
174
175 auto Members() const { return std::tie(yCoeffs, ycbcr); }
176 INLINE_AUTO_MAPPABLE(YuvDesc);
177 };
178
179 struct ColorspaceDesc final {
180 Chromaticities chrom;
181 std::optional<PiecewiseGammaDesc> tf;
182 std::optional<YuvDesc> yuv;
183
184 auto Members() const { return std::tie(chrom, tf, yuv); }
185 INLINE_AUTO_MAPPABLE(ColorspaceDesc);
186 };
187
188 // -
189
190 template <class TT, int NN>
191 struct avec final {
192 using T = TT;
193 static constexpr auto N = NN;
194
195 std::array<T, N> data = {};
196
197 // -
198
199 constexpr avec() = default;
200 constexpr avec(const avec&) = default;
201
202 constexpr avec(const avec<T, N - 1>& v, T a) {
203 for (int i = 0; i < N - 1; i++) {
204 data[i] = v[i];
205 }
206 data[N - 1] = a;
207 }
208 constexpr avec(const avec<T, N - 2>& v, T a, T b) {
209 for (int i = 0; i < N - 2; i++) {
210 data[i] = v[i];
211 }
212 data[N - 2] = a;
213 data[N - 1] = b;
214 }
215
216 MOZ_IMPLICIT constexpr avec(const std::array<T, N>& data) {
217 this->data = data;
218 }
219
220 explicit constexpr avec(const T v) {
221 for (int i = 0; i < N; i++) {
222 data[i] = v;
223 }
224 }
225
226 template <class T2, int N2>
227 explicit constexpr avec(const avec<T2, N2>& v) {
228 const auto n = std::min(N, N2);
229 for (int i = 0; i < n; i++) {
230 data[i] = static_cast<T>(v[i]);
231 }
232 }
233
234 // -
235
236 const auto& operator[](const size_t n) const { return data[n]; }
237 auto& operator[](const size_t n) { return data[n]; }
238
239 template <int i>
240 constexpr auto get() const {
241 return (i < N) ? data[i] : 0;
242 }
243 constexpr auto x() const { return get<0>(); }
244 constexpr auto y() const { return get<1>(); }
245 constexpr auto z() const { return get<2>(); }
246 constexpr auto w() const { return get<3>(); }
247
248 constexpr auto xyz() const { return vec3({x(), y(), z()}); }
249
250 template <int i>
251 void set(const T v) {
252 if (i < N) {
253 data[i] = v;
254 }
255 }
256 void x(const T v) { set<0>(v); }
257 void y(const T v) { set<1>(v); }
258 void z(const T v) { set<2>(v); }
259 void w(const T v) { set<3>(v); }
260
261 // -
262
263 #define _(OP) \
264 friend avec operator OP(const avec a, const avec b) { \
265 avec c; \
266 for (int i = 0; i < N; i++) { \
267 c[i] = a[i] OP b[i]; \
268 } \
269 return c; \
270 } \
271 friend avec operator OP(const avec a, const T b) { \
272 avec c; \
273 for (int i = 0; i < N; i++) { \
274 c[i] = a[i] OP b; \
275 } \
276 return c; \
277 } \
278 friend avec operator OP(const T a, const avec b) { \
279 avec c; \
280 for (int i = 0; i < N; i++) { \
281 c[i] = a OP b[i]; \
282 } \
283 return c; \
284 }
285 _(+)
286 _(-)
287 _(*)
288 _(/)
289 #undef _
290
291 friend bool operator==(const avec a, const avec b) {
292 bool eq = true;
293 for (int i = 0; i < N; i++) {
294 eq &= (a[i] == b[i]);
295 }
296 return eq;
297 }
298 };
299 using vec2 = avec<float, 2>;
300 using vec3 = avec<float, 3>;
301 using vec4 = avec<float, 4>;
302 using ivec3 = avec<int32_t, 3>;
303 using ivec4 = avec<int32_t, 4>;
304
305 template <class T, int N>
306 T dot(const avec<T, N>& a, const avec<T, N>& b) {
307 const auto c = a * b;
308 T ret = 0;
309 for (int i = 0; i < N; i++) {
310 ret += c[i];
311 }
312 return ret;
313 }
314
315 template <class V>
316 V mix(const V& zero, const V& one, const float val) {
317 return zero * (1 - val) + one * val;
318 }
319
320 template <class T, int N>
321 auto min(const avec<T, N>& a, const avec<T, N>& b) {
322 auto ret = avec<T, N>{};
323 for (int i = 0; i < ret.N; i++) {
324 ret[i] = std::min(a[i], b[i]);
325 }
326 return ret;
327 }
328
329 template <class T, int N>
330 auto max(const avec<T, N>& a, const avec<T, N>& b) {
331 auto ret = avec<T, N>{};
332 for (int i = 0; i < ret.N; i++) {
333 ret[i] = std::max(a[i], b[i]);
334 }
335 return ret;
336 }
337
338 template <class T, int N>
339 auto floor(const avec<T, N>& a) {
340 auto ret = avec<T, N>{};
341 for (int i = 0; i < ret.N; i++) {
342 ret[i] = floorf(a[i]);
343 }
344 return ret;
345 }
346
347 template <class T, int N>
348 auto round(const avec<T, N>& a) {
349 auto ret = avec<T, N>{};
350 for (int i = 0; i < ret.N; i++) {
351 ret[i] = roundf(a[i]);
352 }
353 return ret;
354 }
355
356 template <class T, int N>
357 auto abs(const avec<T, N>& a) {
358 auto ret = avec<T, N>{};
359 for (int i = 0; i < ret.N; i++) {
360 ret[i] = std::abs(a[i]);
361 }
362 return ret;
363 }
364
365 // -
366
367 template <int Y_Rows, int X_Cols>
368 struct mat final {
369 static constexpr int y_rows = Y_Rows;
370 static constexpr int x_cols = X_Cols;
371
372 static constexpr auto Identity() {
373 auto ret = mat{};
374 for (int i = 0; i < std::min(x_cols, y_rows); i++) {
375 ret.at(i, i) = 1;
376 }
377 return ret;
378 }
379 static constexpr auto Scale(const avec<float, std::min(x_cols, y_rows)>& v) {
380 auto ret = mat{};
381 for (int i = 0; i < v.N; i++) {
382 ret.at(i, i) = v[i];
383 }
384 return ret;
385 }
386
387 std::array<avec<float, X_Cols>, Y_Rows> rows = {}; // row-major
388
389 // -
390
391 constexpr mat() = default;
392
393 explicit constexpr mat(const std::array<avec<float, X_Cols>, Y_Rows>& rows) {
394 this->rows = rows;
395 }
396
397 template <int Y_Rows2, int X_Cols2>
398 explicit constexpr mat(const mat<Y_Rows2, X_Cols2>& m) {
399 *this = Identity();
400 for (int x = 0; x < std::min(X_Cols, X_Cols2); x++) {
401 for (int y = 0; y < std::min(Y_Rows, Y_Rows2); y++) {
402 at(x, y) = m.at(x, y);
403 }
404 }
405 }
406
407 const auto& at(const int x, const int y) const { return rows.at(y)[x]; }
408 auto& at(const int x, const int y) { return rows.at(y)[x]; }
409
410 friend auto operator*(const mat& a, const avec<float, X_Cols>& b_colvec) {
411 avec<float, Y_Rows> c_colvec;
412 for (int i = 0; i < y_rows; i++) {
413 c_colvec[i] = dot(a.rows.at(i), b_colvec);
414 }
415 return c_colvec;
416 }
417
418 friend auto operator*(const mat& a, const float b) {
419 mat c;
420 for (int x = 0; x < x_cols; x++) {
421 for (int y = 0; y < y_rows; y++) {
422 c.at(x, y) = a.at(x, y) * b;
423 }
424 }
425 return c;
426 }
427 friend auto operator/(const mat& a, const float b) { return a * (1 / b); }
428
429 template <int BCols, int BRows = X_Cols>
430 friend auto operator*(const mat& a, const mat<BRows, BCols>& b) {
431 const auto bt = transpose(b);
432 const auto& b_cols = bt.rows;
433
434 mat<Y_Rows, BCols> c;
435 for (int x = 0; x < BCols; x++) {
436 for (int y = 0; y < Y_Rows; y++) {
437 c.at(x, y) = dot(a.rows.at(y), b_cols.at(x));
438 }
439 }
440 return c;
441 }
442
443 // For e.g. similarity evaluation
444 friend auto operator-(const mat& a, const mat& b) {
445 mat c;
446 for (int y = 0; y < y_rows; y++) {
447 c.rows[y] = a.rows[y] - b.rows[y];
448 }
449 return c;
450 }
451 };
452
453 template <class M>
454 inline float dotDifference(const M& a, const M& b) {
455 const auto c = a - b;
456 const auto d = c * avec<float, M::x_cols>(1);
457 const auto d2 = dot(d, d);
458 return d2;
459 }
460 template <class M>
461 inline bool approx(const M& a, const M& b, const float eps = 0.0001) {
462 const auto errSquared = dotDifference(a, b);
463 return errSquared <= (eps * eps);
464 }
465
466 using mat3 = mat<3, 3>;
467 using mat4 = mat<4, 4>;
468
469 inline float determinant(const mat<1, 1>& m) { return m.at(0, 0); }
470 template <class T>
471 float determinant(const T& m) {
472 static_assert(T::x_cols == T::y_rows);
473
474 float ret = 0;
475 for (int i = 0; i < T::x_cols; i++) {
476 const auto cofact = cofactor(m, i, 0);
477 ret += m.at(i, 0) * cofact;
478 }
479 return ret;
480 }
481
482 // -
483
484 template <class T>
485 float cofactor(const T& m, const int x_col, const int y_row) {
486 ASSERT(0 <= x_col && x_col < T::x_cols);
487 ASSERT(0 <= y_row && y_row < T::y_rows);
488
489 auto cofactor = minor_val(m, x_col, y_row);
490 if ((x_col + y_row) % 2 == 1) {
491 cofactor *= -1;
492 }
493 return cofactor;
494 }
495
496 // -
497
498 // Unfortunately, can't call this `minor(...)` because there is
499 // `#define minor(dev) gnu_dev_minor (dev)`
500 // in /usr/include/x86_64-linux-gnu/sys/sysmacros.h:62
501 template <class T>
502 float minor_val(const T& a, const int skip_x, const int skip_y) {
503 ASSERT(0 <= skip_x && skip_x < T::x_cols);
504 ASSERT(0 <= skip_y && skip_y < T::y_rows);
505
506 // A minor matrix is a matrix without its x_col and y_row.
507 mat<T::y_rows - 1, T::x_cols - 1> b;
508
509 int x_skips = 0;
510 for (int ax = 0; ax < T::x_cols; ax++) {
511 if (ax == skip_x) {
512 x_skips = 1;
513 continue;
514 }
515
516 int y_skips = 0;
517 for (int ay = 0; ay < T::y_rows; ay++) {
518 if (ay == skip_y) {
519 y_skips = 1;
520 continue;
521 }
522
523 b.at(ax - x_skips, ay - y_skips) = a.at(ax, ay);
524 }
525 }
526
527 const auto minor = determinant(b);
528 return minor;
529 }
530
531 // -
532
533 /// The matrix of cofactors.
534 template <class T>
535 auto comatrix(const T& a) {
536 auto b = T{};
537 for (int x = 0; x < T::x_cols; x++) {
538 for (int y = 0; y < T::y_rows; y++) {
539 b.at(x, y) = cofactor(a, x, y);
540 }
541 }
542 return b;
543 }
544
545 // -
546
547 template <class T>
548 auto transpose(const T& a) {
549 auto b = mat<T::x_cols, T::y_rows>{};
550 for (int x = 0; x < T::x_cols; x++) {
551 for (int y = 0; y < T::y_rows; y++) {
552 b.at(y, x) = a.at(x, y);
553 }
554 }
555 return b;
556 }
557
558 // -
559
560 template <class T>
561 inline T inverse(const T& a) {
562 const auto det = determinant(a);
563 const auto comat = comatrix(a);
564 const auto adjugate = transpose(comat);
565 const auto inv = adjugate / det;
566 return inv;
567 }
568
569 // -
570
571 template <class F>
572 void ForEachIntWithin(const ivec3 size, const F& f) {
573 ivec3 p;
574 for (p.z(0); p.z() < size.z(); p.z(p.z() + 1)) {
575 for (p.y(0); p.y() < size.y(); p.y(p.y() + 1)) {
576 for (p.x(0); p.x() < size.x(); p.x(p.x() + 1)) {
577 f(p);
578 }
579 }
580 }
581 }
582 template <class F>
583 void ForEachSampleWithin(const ivec3 size, const F& f) {
584 const auto div = vec3(size - 1);
585 ForEachIntWithin(size, [&](const ivec3& isrc) {
586 const auto fsrc = vec3(isrc) / div;
587 f(fsrc);
588 });
589 }
590
591 // -
592
593 struct Lut3 final {
594 ivec3 size;
595 std::vector<vec3> data;
596
597 // -
598
599 static Lut3 Create(const ivec3 size) {
600 Lut3 lut;
601 lut.size = size;
602 lut.data.resize(size.x() * size.y() * size.z());
603 return lut;
604 }
605
606 // -
607
608 /// p: [0, N-1] (clamps)
609 size_t Index(ivec3 p) const {
610 const auto scales = ivec3({1, size.x(), size.x() * size.y()});
611 p = max(ivec3(0), min(p, size - 1)); // clamp
612 return dot(p, scales);
613 }
614
615 // -
616
617 template <class F>
618 void SetMap(const F& dstFromSrc01) {
619 ForEachIntWithin(size, [&](const ivec3 p) {
620 const auto i = Index(p);
621 const auto src01 = vec3(p) / vec3(size - 1);
622 const auto dstVal = dstFromSrc01(src01);
623 data.at(i) = dstVal;
624 });
625 }
626
627 // -
628
629 /// p: [0, N-1] (clamps)
630 vec3 Fetch(ivec3 p) const {
631 const auto i = Index(p);
632 return data.at(i);
633 }
634
635 /// in01: [0.0, 1.0] (clamps)
636 vec3 Sample(vec3 in01) const;
637 };
638
639 // -
640
641 /**
642 Naively, it would be ideal to map directly from ycbcr to rgb,
643 but headroom and footroom are problematic: For e.g. narrow-range-8-bit,
644 our naive LUT would start at absolute y=0/255. However, values only start
645 at y=16/255, and depending on where your first LUT sample is, you might get
646 very poor approximations for y=16/255.
647 Further, even for full-range-8-bit, y=-0.5 is encoded as 1/255. U and v
648 aren't *as* important as y, but we should try be accurate for the min and
649 max values. Additionally, it would be embarassing to get whites/greys wrong,
650 so preserving u=0.0 should also be a goal.
651 Finally, when using non-linear transfer functions, the linear approximation of a
652 point between two samples will be fairly inaccurate.
653 We preserve min and max by choosing our input range such that min and max are
654 the endpoints of their LUT axis.
655 We preserve accuracy (at and around) mid by choosing odd sizes for dimentions.
656
657 But also, the LUT is surprisingly robust, so check if the simple version works
658 before adding complexity!
659 **/
660
661 struct ColorspaceTransform final {
662 ColorspaceDesc srcSpace;
663 ColorspaceDesc dstSpace;
664 mat4 srcRgbTfFromSrc;
665 std::optional<PiecewiseGammaDesc> srcTf;
666 mat3 dstRgbLinFromSrcRgbLin;
667 std::optional<PiecewiseGammaDesc> dstTf;
668 mat4 dstFromDstRgbTf;
669
670 static ColorspaceTransform Create(const ColorspaceDesc& src,
671 const ColorspaceDesc& dst);
672
673 // -
674
675 vec3 DstFromSrc(vec3 src) const;
676
677 std::optional<mat4> ToMat4() const;
678
679 Lut3 ToLut3(const ivec3 size) const;
680 Lut3 ToLut3() const {
681 auto defaultSize = ivec3({31, 31, 15}); // Order of importance: G, R, B
682 if (srcSpace.yuv) {
683 defaultSize = ivec3({31, 15, 31}); // Y, Cb, Cr
684 }
685 return ToLut3(defaultSize);
686 }
687 };
688
689 // -
690
691 struct RgbTransferTables {
692 std::vector<float> r;
693 std::vector<float> g;
694 std::vector<float> b;
695 };
696 float GuessGamma(const std::vector<float>& vals, float exp_guess = 1.0);
697
698 static constexpr auto D65 = vec2{{0.3127, 0.3290}};
699 static constexpr auto D50 = vec2{{0.34567, 0.35850}};
700 mat3 XyzAFromXyzB_BradfordLinear(const vec2 xyA, const vec2 xyB);
701
702 // -
703
704 struct ColorProfileDesc {
705 // ICC profiles are phrased as PCS-from-encoded (PCS is CIEXYZ-D50)
706 // However, all of our colorspaces are D65, so let's normalize to that,
707 // even though it's a reversible transform.
708 color::mat4 rgbFromYcbcr = color::mat4::Identity();
709 RgbTransferTables linearFromTf;
710 color::mat3 xyzd65FromLinearRgb = color::mat3::Identity();
711
712 static ColorProfileDesc From(const ColorspaceDesc&);
713 static ColorProfileDesc From(const qcms_profile&);
714 };
715
716 template <class C>
717 inline float SampleOutByIn(const C& outByIn, const float in) {
718 switch (outByIn.size()) {
719 case 0:
720 return in;
721 case 1:
722 return outByIn.at(0);
723 }
724 MOZ_ASSERT(outByIn.size() >= 2);
725 const auto begin = outByIn.begin();
726
727 const auto in0i = size_t(floorf(in * (outByIn.size() - 1)));
728 const auto out0_itr = begin + std::min(in0i, outByIn.size() - 2);
729
730 const auto in0 = float(out0_itr - begin) / (outByIn.size() - 1);
731 const auto out0 = *out0_itr;
732 const auto d_in = float(1) / (outByIn.size() - 1);
733 const auto d_out = *(out0_itr + 1) - *out0_itr;
734
735 const auto out = out0 + (d_out / d_in) * (in - in0);
736 // printf("SampleOutByIn(%f)->%f\n", in, out);
737 return out;
738 }
739
740 template <class C>
741 inline float SampleInByOut(const C& outByIn, const float out) {
742 MOZ_ASSERT(outByIn.size() >= 2);
743 const auto begin = outByIn.begin();
744
745 const auto out0_itr = std::lower_bound(begin + 1, outByIn.end() - 1, out) - 1;
746
747 const auto in0 = float(out0_itr - begin) / (outByIn.size() - 1);
748 const auto out0 = *out0_itr;
749 const auto d_in = float(1) / (outByIn.size() - 1);
750 const auto d_out = *(out0_itr + 1) - *out0_itr;
751
752 // printf("%f + (%f / %f) * (%f - %f)\n", in0, d_in, d_out, out, out0);
753 const auto in = in0 + (d_in / d_out) * (out - out0);
754 // printf("SampleInByOut(%f)->%f\n", out, in);
755 return in;
756 }
757
758 template <class C, class FnLessEqualT = std::less_equal<typename C::value_type>>
759 inline bool IsMonotonic(const C& vals, const FnLessEqualT& LessEqual = {}) {
760 bool ok = true;
761 const auto begin = vals.begin();
762 for (size_t i = 1; i < vals.size(); i++) {
763 const auto itr = begin + i;
764 ok &= LessEqual(*(itr - 1), *itr);
765 // Assert(true, [&]() {
766 // return prints("[%zu]->%f <= [%zu]->%f", i-1, *(itr-1), i, *itr);
767 // });
768 }
769 return ok;
770 }
771
772 template <class T, class I>
773 inline std::optional<I> SeekNeq(const T& ref, const I first, const I last) {
774 const auto inc = (last - first) > 0 ? 1 : -1;
775 auto itr = first;
776 while (true) {
777 if (*itr != ref) return itr;
778 if (itr == last) return {};
779 itr += inc;
780 }
781 }
782
783 template <class T>
784 struct TwoPoints {
785 struct {
786 T x;
787 T y;
788 } p0;
789 struct {
790 T x;
791 T y;
792 } p1;
793
794 T y(const T x) const {
795 const auto dx = p1.x - p0.x;
796 const auto dy = p1.y - p0.y;
797 return p0.y + dy / dx * (x - p0.x);
798 }
799 };
800
801 /// Fills `vals` with `x:[0..vals.size()-1] => line.y(x)`.
802 template <class T>
803 static void LinearFill(T& vals, const TwoPoints<float>& line) {
804 float x = -1;
805 for (auto& val : vals) {
806 x += 1;
807 val = line.y(x);
808 }
809 }
810
811 // -
812
813 inline void DequantizeMonotonic(const Span<float> vals) {
814 MOZ_ASSERT(IsMonotonic(vals));
815
816 const auto first = vals.begin();
817 const auto end = vals.end();
818 if (first == end) return;
819 const auto last = end - 1;
820 if (first == last) return;
821
822 // Three monotonic cases:
823 // 1. [0,0,0,0]
824 // 2. [0,0,1,1]
825 // 3. [0,1,1,2]
826
827 const auto body_first = SeekNeq(*first, first, last);
828 if (!body_first) {
829 // E.g. [0,0,0,0]
830 return;
831 }
832
833 const auto body_last = SeekNeq(*last, last, *body_first);
834 if (!body_last) {
835 // E.g. [0,0,1,1]
836 // This isn't the most accurate, but close enough.
837 // print("#2: %s", to_str(vals).c_str());
838 LinearFill(vals, {
839 {0, *first},
840 {float(vals.size() - 1), *last},
841 });
842 // print(" -> %s\n", to_str(vals).c_str());
843 return;
844 }
845
846 // E.g. [0,1,1,2]
847 // ^^^ body
848 // => f(0.5)->0.5, f(2.5)->1.5
849 // => f(x) = f(x0) + (x-x0) * (f(x1) - f(x0)) / (x1-x0)
850 // => f(x) = f(x0) + (x-x0) * dfdx
851
852 const auto head_end = *body_first;
853 const auto head = vals.subspan(0, head_end - vals.begin());
854 const auto tail_begin = *body_last + 1;
855 const auto tail = vals.subspan(tail_begin - vals.begin());
856 // print("head tail: %s %s\n",
857 // to_str(head).c_str(),
858 // to_str(tail).c_str());
859
860 // const auto body = vals->subspan(head.size(), vals->size()-tail.size());
861 auto next_part_first = head_end;
862 while (next_part_first != tail_begin) {
863 const auto part_first = next_part_first;
864 // print("part_first: %f\n", *part_first);
865 next_part_first = *SeekNeq(*part_first, part_first, tail_begin);
866 // print("next_part_first: %f\n", *next_part_first);
867 const auto part =
868 Span<float>{part_first, size_t(next_part_first - part_first)};
869 // print("part: %s\n", to_str(part).c_str());
870 const auto prev_part_last = part_first - 1;
871 const auto part_last = next_part_first - 1;
872 const auto line = TwoPoints<float>{
873 {-0.5, (*prev_part_last + *part_first) / 2},
874 {part.size() - 0.5f, (*part_last + *next_part_first) / 2},
875 };
876 LinearFill(part, line);
877 }
878
879 static constexpr bool INFER_HEAD_TAIL_FROM_BODY_EDGE = false;
880 // Basically ignore contents of head and tail, and infer from edges of body.
881 // print("3: %s\n", to_str(vals).c_str());
882 if (!IsMonotonic(head, std::less<float>{})) {
883 if (!INFER_HEAD_TAIL_FROM_BODY_EDGE) {
884 LinearFill(head,
885 {
886 {0, *head.begin()},
887 {head.size() - 0.5f, (*(head.end() - 1) + *head_end) / 2},
888 });
889 } else {
890 LinearFill(head, {
891 {head.size() + 0.0f, *head_end},
892 {head.size() + 1.0f, *(head_end + 1)},
893 });
894 }
895 }
896 if (!IsMonotonic(tail, std::less<float>{})) {
897 if (!INFER_HEAD_TAIL_FROM_BODY_EDGE) {
898 LinearFill(tail, {
899 {-0.5, (*(tail_begin - 1) + *tail.begin()) / 2},
900 {tail.size() - 1.0f, *(tail.end() - 1)},
901 });
902 } else {
903 LinearFill(tail, {
904 {-2.0f, *(tail_begin - 2)},
905 {-1.0f, *(tail_begin - 1)},
906 });
907 }
908 }
909 // print("3: %s\n", to_str(vals).c_str());
910 MOZ_ASSERT(IsMonotonic(vals, std::less<float>{}));
911
912 // Rescale, because we tend to lose range.
913 static constexpr bool RESCALE = false;
914 if (RESCALE) {
915 const auto firstv = *first;
916 const auto lastv = *last;
917 for (auto& val : vals) {
918 val = (val - firstv) / (lastv - firstv);
919 }
920 }
921 // print("4: %s\n", to_str(vals).c_str());
922 }
923
924 template <class In, class Out>
925 static void InvertLut(const In& lut, Out* const out_invertedLut) {
926 MOZ_ASSERT(IsMonotonic(lut));
927 auto plut = &lut;
928 auto vec = std::vector<float>{};
929 if (!IsMonotonic(lut, std::less<float>{})) {
930 // print("Not strictly monotonic...\n");
931 vec.assign(lut.begin(), lut.end());
932 DequantizeMonotonic(vec);
933 plut = &vec;
934 // print(" Now strictly monotonic: %i: %s\n",
935 // int(IsMonotonic(*plut, std::less<float>{})), to_str(*plut).c_str());
936 MOZ_ASSERT(IsMonotonic(*plut, std::less<float>{}));
937 }
938 MOZ_ASSERT(plut->size() >= 2);
939
940 auto& ret = *out_invertedLut;
941 for (size_t i_out = 0; i_out < ret.size(); i_out++) {
942 const auto f_out = i_out / float(ret.size() - 1);
943 const auto f_in = SampleInByOut(*plut, f_out);
944 ret[i_out] = f_in;
945 }
946
947 MOZ_ASSERT(IsMonotonic(ret));
948 MOZ_ASSERT(IsMonotonic(ret, std::less<float>{}));
949 }
950
951 // -
952
953 struct ColorProfileConversionDesc {
954 // ICC profiles are phrased as PCS-from-encoded (PCS is CIEXYZ-D50)
955 color::mat4 srcRgbFromSrcYuv = color::mat4::Identity();
956 RgbTransferTables srcLinearFromSrcTf;
957 color::mat3 dstLinearFromSrcLinear = color::mat3::Identity();
958 RgbTransferTables dstTfFromDstLinear;
959
960 struct FromDesc {
961 ColorProfileDesc src;
962 ColorProfileDesc dst;
963 };
964 static ColorProfileConversionDesc From(const FromDesc&);
965
966 vec3 Apply(const vec3 src) const {
967 const auto srcRgb = vec3(srcRgbFromSrcYuv * vec4(src, 1));
968 const auto srcLinear = vec3{{
969 SampleOutByIn(srcLinearFromSrcTf.r, srcRgb.x()),
970 SampleOutByIn(srcLinearFromSrcTf.g, srcRgb.y()),
971 SampleOutByIn(srcLinearFromSrcTf.b, srcRgb.z()),
972 }};
973 const auto dstLinear = dstLinearFromSrcLinear * srcLinear;
974 const auto dstRgb = vec3{{
975 SampleOutByIn(dstTfFromDstLinear.r, dstLinear.x()),
976 SampleOutByIn(dstTfFromDstLinear.g, dstLinear.y()),
977 SampleOutByIn(dstTfFromDstLinear.b, dstLinear.z()),
978 }};
979 return dstRgb;
980 }
981 };
982
983 } // namespace mozilla::color
984
985 #undef ASSERT
986
987 #endif // MOZILLA_GFX_GL_COLORSPACES_H_