summaryrefslogtreecommitdiffstats
path: root/gfx/gl/Colorspaces.h
diff options
context:
space:
mode:
Diffstat (limited to 'gfx/gl/Colorspaces.h')
-rw-r--r--gfx/gl/Colorspaces.h987
1 files changed, 987 insertions, 0 deletions
diff --git a/gfx/gl/Colorspaces.h b/gfx/gl/Colorspaces.h
new file mode 100644
index 0000000000..8f36854d2d
--- /dev/null
+++ b/gfx/gl/Colorspaces.h
@@ -0,0 +1,987 @@
+/* This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
+
+#ifndef MOZILLA_GFX_GL_COLORSPACES_H_
+#define MOZILLA_GFX_GL_COLORSPACES_H_
+
+// Reference: https://hackmd.io/0wkiLmP7RWOFjcD13M870A
+
+// We are going to be doing so, so many transforms, so descriptive labels are
+// critical.
+
+// Colorspace background info: https://hackmd.io/0wkiLmP7RWOFjcD13M870A
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <cstdint>
+#include <cstdlib>
+#include <functional>
+#include <optional>
+#include <vector>
+
+#include "AutoMappable.h"
+#include "mozilla/Assertions.h"
+#include "mozilla/Attributes.h"
+#include "mozilla/Span.h"
+
+#ifdef DEBUG
+# define ASSERT(EXPR) \
+ do { \
+ if (!(EXPR)) { \
+ __builtin_trap(); \
+ } \
+ } while (false)
+#else
+# define ASSERT(EXPR) (void)(EXPR)
+#endif
+
+struct _qcms_profile;
+typedef struct _qcms_profile qcms_profile;
+
+namespace mozilla::color {
+
+struct YuvLumaCoeffs final {
+ float r = 0.2126;
+ float g = 0.7152;
+ float b = 0.0722;
+
+ auto Members() const { return std::tie(r, g, b); }
+ INLINE_AUTO_MAPPABLE(YuvLumaCoeffs)
+
+ static constexpr auto Rec709() { return YuvLumaCoeffs(); }
+
+ static constexpr auto Rec2020() {
+ return YuvLumaCoeffs{0.2627, 0.6780, 0.0593};
+ }
+};
+
+struct PiecewiseGammaDesc final {
+ // tf = { k * linear | linear < b
+ // { a * pow(linear, 1/g) - (1-a) | linear >= b
+
+ // Default to Srgb
+ float a = 1.055;
+ float b = 0.04045 / 12.92;
+ float g = 2.4;
+ float k = 12.92;
+
+ auto Members() const { return std::tie(a, b, g, k); }
+ INLINE_AUTO_MAPPABLE(PiecewiseGammaDesc)
+
+ static constexpr auto Srgb() { return PiecewiseGammaDesc(); }
+ static constexpr auto DisplayP3() { return Srgb(); }
+
+ static constexpr auto Rec709() {
+ return PiecewiseGammaDesc{
+ 1.099,
+ 0.018,
+ 1.0 / 0.45, // ~2.222
+ 4.5,
+ };
+ }
+ // FYI: static constexpr auto Rec2020_10bit() { return Rec709(); }
+ static constexpr auto Rec2020_12bit() {
+ return PiecewiseGammaDesc{
+ 1.0993,
+ 0.0181,
+ 1.0 / 0.45, // ~2.222
+ 4.5,
+ };
+ }
+};
+
+struct YcbcrDesc final {
+ float y0 = 16 / 255.0;
+ float y1 = 235 / 255.0;
+ float u0 = 128 / 255.0;
+ float uPlusHalf = 240 / 255.0;
+
+ auto Members() const { return std::tie(y0, y1, u0, uPlusHalf); }
+ INLINE_AUTO_MAPPABLE(YcbcrDesc)
+
+ static constexpr auto Narrow8() { // AKA limited/studio/tv
+ return YcbcrDesc();
+ }
+ static constexpr auto Full8() { // AKA pc
+ return YcbcrDesc{
+ 0 / 255.0,
+ 255 / 255.0,
+ 128 / 255.0,
+ 254 / 255.0,
+ };
+ }
+ static constexpr auto Float() { // Best for a LUT
+ return YcbcrDesc{0.0, 1.0, 0.5, 1.0};
+ }
+};
+
+struct Chromaticities final {
+ float rx = 0.640;
+ float ry = 0.330;
+ float gx = 0.300;
+ float gy = 0.600;
+ float bx = 0.150;
+ float by = 0.060;
+ // D65:
+ static constexpr float wx = 0.3127;
+ static constexpr float wy = 0.3290;
+
+ auto Members() const { return std::tie(rx, ry, gx, gy, bx, by); }
+ INLINE_AUTO_MAPPABLE(Chromaticities)
+
+ // -
+
+ static constexpr auto Rec709() { // AKA limited/studio/tv
+ return Chromaticities();
+ }
+ static constexpr auto Srgb() { return Rec709(); }
+
+ static constexpr auto Rec601_625_Pal() {
+ auto ret = Rec709();
+ ret.gx = 0.290;
+ return ret;
+ }
+ static constexpr auto Rec601_525_Ntsc() {
+ return Chromaticities{
+ 0.630, 0.340, // r
+ 0.310, 0.595, // g
+ 0.155, 0.070, // b
+ };
+ }
+ static constexpr auto Rec2020() {
+ return Chromaticities{
+ 0.708, 0.292, // r
+ 0.170, 0.797, // g
+ 0.131, 0.046, // b
+ };
+ }
+ static constexpr auto DisplayP3() {
+ return Chromaticities{
+ 0.680, 0.320, // r
+ 0.265, 0.690, // g
+ 0.150, 0.060, // b
+ };
+ }
+};
+
+// -
+
+struct YuvDesc final {
+ YuvLumaCoeffs yCoeffs;
+ YcbcrDesc ycbcr;
+
+ auto Members() const { return std::tie(yCoeffs, ycbcr); }
+ INLINE_AUTO_MAPPABLE(YuvDesc);
+};
+
+struct ColorspaceDesc final {
+ Chromaticities chrom;
+ std::optional<PiecewiseGammaDesc> tf;
+ std::optional<YuvDesc> yuv;
+
+ auto Members() const { return std::tie(chrom, tf, yuv); }
+ INLINE_AUTO_MAPPABLE(ColorspaceDesc);
+};
+
+// -
+
+template <class TT, int NN>
+struct avec final {
+ using T = TT;
+ static constexpr auto N = NN;
+
+ std::array<T, N> data = {};
+
+ // -
+
+ constexpr avec() = default;
+ constexpr avec(const avec&) = default;
+
+ constexpr avec(const avec<T, N - 1>& v, T a) {
+ for (int i = 0; i < N - 1; i++) {
+ data[i] = v[i];
+ }
+ data[N - 1] = a;
+ }
+ constexpr avec(const avec<T, N - 2>& v, T a, T b) {
+ for (int i = 0; i < N - 2; i++) {
+ data[i] = v[i];
+ }
+ data[N - 2] = a;
+ data[N - 1] = b;
+ }
+
+ MOZ_IMPLICIT constexpr avec(const std::array<T, N>& data) {
+ this->data = data;
+ }
+
+ explicit constexpr avec(const T v) {
+ for (int i = 0; i < N; i++) {
+ data[i] = v;
+ }
+ }
+
+ template <class T2, int N2>
+ explicit constexpr avec(const avec<T2, N2>& v) {
+ const auto n = std::min(N, N2);
+ for (int i = 0; i < n; i++) {
+ data[i] = static_cast<T>(v[i]);
+ }
+ }
+
+ // -
+
+ const auto& operator[](const size_t n) const { return data[n]; }
+ auto& operator[](const size_t n) { return data[n]; }
+
+ template <int i>
+ constexpr auto get() const {
+ return (i < N) ? data[i] : 0;
+ }
+ constexpr auto x() const { return get<0>(); }
+ constexpr auto y() const { return get<1>(); }
+ constexpr auto z() const { return get<2>(); }
+ constexpr auto w() const { return get<3>(); }
+
+ constexpr auto xyz() const { return vec3({x(), y(), z()}); }
+
+ template <int i>
+ void set(const T v) {
+ if (i < N) {
+ data[i] = v;
+ }
+ }
+ void x(const T v) { set<0>(v); }
+ void y(const T v) { set<1>(v); }
+ void z(const T v) { set<2>(v); }
+ void w(const T v) { set<3>(v); }
+
+ // -
+
+#define _(OP) \
+ friend avec operator OP(const avec a, const avec b) { \
+ avec c; \
+ for (int i = 0; i < N; i++) { \
+ c[i] = a[i] OP b[i]; \
+ } \
+ return c; \
+ } \
+ friend avec operator OP(const avec a, const T b) { \
+ avec c; \
+ for (int i = 0; i < N; i++) { \
+ c[i] = a[i] OP b; \
+ } \
+ return c; \
+ } \
+ friend avec operator OP(const T a, const avec b) { \
+ avec c; \
+ for (int i = 0; i < N; i++) { \
+ c[i] = a OP b[i]; \
+ } \
+ return c; \
+ }
+ _(+)
+ _(-)
+ _(*)
+ _(/)
+#undef _
+
+ friend bool operator==(const avec a, const avec b) {
+ bool eq = true;
+ for (int i = 0; i < N; i++) {
+ eq &= (a[i] == b[i]);
+ }
+ return eq;
+ }
+};
+using vec2 = avec<float, 2>;
+using vec3 = avec<float, 3>;
+using vec4 = avec<float, 4>;
+using ivec3 = avec<int32_t, 3>;
+using ivec4 = avec<int32_t, 4>;
+
+template <class T, int N>
+T dot(const avec<T, N>& a, const avec<T, N>& b) {
+ const auto c = a * b;
+ T ret = 0;
+ for (int i = 0; i < N; i++) {
+ ret += c[i];
+ }
+ return ret;
+}
+
+template <class V>
+V mix(const V& zero, const V& one, const float val) {
+ return zero * (1 - val) + one * val;
+}
+
+template <class T, int N>
+auto min(const avec<T, N>& a, const avec<T, N>& b) {
+ auto ret = avec<T, N>{};
+ for (int i = 0; i < ret.N; i++) {
+ ret[i] = std::min(a[i], b[i]);
+ }
+ return ret;
+}
+
+template <class T, int N>
+auto max(const avec<T, N>& a, const avec<T, N>& b) {
+ auto ret = avec<T, N>{};
+ for (int i = 0; i < ret.N; i++) {
+ ret[i] = std::max(a[i], b[i]);
+ }
+ return ret;
+}
+
+template <class T, int N>
+auto floor(const avec<T, N>& a) {
+ auto ret = avec<T, N>{};
+ for (int i = 0; i < ret.N; i++) {
+ ret[i] = floorf(a[i]);
+ }
+ return ret;
+}
+
+template <class T, int N>
+auto round(const avec<T, N>& a) {
+ auto ret = avec<T, N>{};
+ for (int i = 0; i < ret.N; i++) {
+ ret[i] = roundf(a[i]);
+ }
+ return ret;
+}
+
+template <class T, int N>
+auto abs(const avec<T, N>& a) {
+ auto ret = avec<T, N>{};
+ for (int i = 0; i < ret.N; i++) {
+ ret[i] = std::abs(a[i]);
+ }
+ return ret;
+}
+
+// -
+
+template <int Y_Rows, int X_Cols>
+struct mat final {
+ static constexpr int y_rows = Y_Rows;
+ static constexpr int x_cols = X_Cols;
+
+ static constexpr auto Identity() {
+ auto ret = mat{};
+ for (int i = 0; i < std::min(x_cols, y_rows); i++) {
+ ret.at(i, i) = 1;
+ }
+ return ret;
+ }
+ static constexpr auto Scale(const avec<float, std::min(x_cols, y_rows)>& v) {
+ auto ret = mat{};
+ for (int i = 0; i < v.N; i++) {
+ ret.at(i, i) = v[i];
+ }
+ return ret;
+ }
+
+ std::array<avec<float, X_Cols>, Y_Rows> rows = {}; // row-major
+
+ // -
+
+ constexpr mat() = default;
+
+ explicit constexpr mat(const std::array<avec<float, X_Cols>, Y_Rows>& rows) {
+ this->rows = rows;
+ }
+
+ template <int Y_Rows2, int X_Cols2>
+ explicit constexpr mat(const mat<Y_Rows2, X_Cols2>& m) {
+ *this = Identity();
+ for (int x = 0; x < std::min(X_Cols, X_Cols2); x++) {
+ for (int y = 0; y < std::min(Y_Rows, Y_Rows2); y++) {
+ at(x, y) = m.at(x, y);
+ }
+ }
+ }
+
+ const auto& at(const int x, const int y) const { return rows.at(y)[x]; }
+ auto& at(const int x, const int y) { return rows.at(y)[x]; }
+
+ friend auto operator*(const mat& a, const avec<float, X_Cols>& b_colvec) {
+ avec<float, Y_Rows> c_colvec;
+ for (int i = 0; i < y_rows; i++) {
+ c_colvec[i] = dot(a.rows.at(i), b_colvec);
+ }
+ return c_colvec;
+ }
+
+ friend auto operator*(const mat& a, const float b) {
+ mat c;
+ for (int x = 0; x < x_cols; x++) {
+ for (int y = 0; y < y_rows; y++) {
+ c.at(x, y) = a.at(x, y) * b;
+ }
+ }
+ return c;
+ }
+ friend auto operator/(const mat& a, const float b) { return a * (1 / b); }
+
+ template <int BCols, int BRows = X_Cols>
+ friend auto operator*(const mat& a, const mat<BRows, BCols>& b) {
+ const auto bt = transpose(b);
+ const auto& b_cols = bt.rows;
+
+ mat<Y_Rows, BCols> c;
+ for (int x = 0; x < BCols; x++) {
+ for (int y = 0; y < Y_Rows; y++) {
+ c.at(x, y) = dot(a.rows.at(y), b_cols.at(x));
+ }
+ }
+ return c;
+ }
+
+ // For e.g. similarity evaluation
+ friend auto operator-(const mat& a, const mat& b) {
+ mat c;
+ for (int y = 0; y < y_rows; y++) {
+ c.rows[y] = a.rows[y] - b.rows[y];
+ }
+ return c;
+ }
+};
+
+template <class M>
+inline float dotDifference(const M& a, const M& b) {
+ const auto c = a - b;
+ const auto d = c * avec<float, M::x_cols>(1);
+ const auto d2 = dot(d, d);
+ return d2;
+}
+template <class M>
+inline bool approx(const M& a, const M& b, const float eps = 0.0001) {
+ const auto errSquared = dotDifference(a, b);
+ return errSquared <= (eps * eps);
+}
+
+using mat3 = mat<3, 3>;
+using mat4 = mat<4, 4>;
+
+inline float determinant(const mat<1, 1>& m) { return m.at(0, 0); }
+template <class T>
+float determinant(const T& m) {
+ static_assert(T::x_cols == T::y_rows);
+
+ float ret = 0;
+ for (int i = 0; i < T::x_cols; i++) {
+ const auto cofact = cofactor(m, i, 0);
+ ret += m.at(i, 0) * cofact;
+ }
+ return ret;
+}
+
+// -
+
+template <class T>
+float cofactor(const T& m, const int x_col, const int y_row) {
+ ASSERT(0 <= x_col && x_col < T::x_cols);
+ ASSERT(0 <= y_row && y_row < T::y_rows);
+
+ auto cofactor = minor_val(m, x_col, y_row);
+ if ((x_col + y_row) % 2 == 1) {
+ cofactor *= -1;
+ }
+ return cofactor;
+}
+
+// -
+
+// Unfortunately, can't call this `minor(...)` because there is
+// `#define minor(dev) gnu_dev_minor (dev)`
+// in /usr/include/x86_64-linux-gnu/sys/sysmacros.h:62
+template <class T>
+float minor_val(const T& a, const int skip_x, const int skip_y) {
+ ASSERT(0 <= skip_x && skip_x < T::x_cols);
+ ASSERT(0 <= skip_y && skip_y < T::y_rows);
+
+ // A minor matrix is a matrix without its x_col and y_row.
+ mat<T::y_rows - 1, T::x_cols - 1> b;
+
+ int x_skips = 0;
+ for (int ax = 0; ax < T::x_cols; ax++) {
+ if (ax == skip_x) {
+ x_skips = 1;
+ continue;
+ }
+
+ int y_skips = 0;
+ for (int ay = 0; ay < T::y_rows; ay++) {
+ if (ay == skip_y) {
+ y_skips = 1;
+ continue;
+ }
+
+ b.at(ax - x_skips, ay - y_skips) = a.at(ax, ay);
+ }
+ }
+
+ const auto minor = determinant(b);
+ return minor;
+}
+
+// -
+
+/// The matrix of cofactors.
+template <class T>
+auto comatrix(const T& a) {
+ auto b = T{};
+ for (int x = 0; x < T::x_cols; x++) {
+ for (int y = 0; y < T::y_rows; y++) {
+ b.at(x, y) = cofactor(a, x, y);
+ }
+ }
+ return b;
+}
+
+// -
+
+template <class T>
+auto transpose(const T& a) {
+ auto b = mat<T::x_cols, T::y_rows>{};
+ for (int x = 0; x < T::x_cols; x++) {
+ for (int y = 0; y < T::y_rows; y++) {
+ b.at(y, x) = a.at(x, y);
+ }
+ }
+ return b;
+}
+
+// -
+
+template <class T>
+inline T inverse(const T& a) {
+ const auto det = determinant(a);
+ const auto comat = comatrix(a);
+ const auto adjugate = transpose(comat);
+ const auto inv = adjugate / det;
+ return inv;
+}
+
+// -
+
+template <class F>
+void ForEachIntWithin(const ivec3 size, const F& f) {
+ ivec3 p;
+ for (p.z(0); p.z() < size.z(); p.z(p.z() + 1)) {
+ for (p.y(0); p.y() < size.y(); p.y(p.y() + 1)) {
+ for (p.x(0); p.x() < size.x(); p.x(p.x() + 1)) {
+ f(p);
+ }
+ }
+ }
+}
+template <class F>
+void ForEachSampleWithin(const ivec3 size, const F& f) {
+ const auto div = vec3(size - 1);
+ ForEachIntWithin(size, [&](const ivec3& isrc) {
+ const auto fsrc = vec3(isrc) / div;
+ f(fsrc);
+ });
+}
+
+// -
+
+struct Lut3 final {
+ ivec3 size;
+ std::vector<vec3> data;
+
+ // -
+
+ static Lut3 Create(const ivec3 size) {
+ Lut3 lut;
+ lut.size = size;
+ lut.data.resize(size.x() * size.y() * size.z());
+ return lut;
+ }
+
+ // -
+
+ /// p: [0, N-1] (clamps)
+ size_t Index(ivec3 p) const {
+ const auto scales = ivec3({1, size.x(), size.x() * size.y()});
+ p = max(ivec3(0), min(p, size - 1)); // clamp
+ return dot(p, scales);
+ }
+
+ // -
+
+ template <class F>
+ void SetMap(const F& dstFromSrc01) {
+ ForEachIntWithin(size, [&](const ivec3 p) {
+ const auto i = Index(p);
+ const auto src01 = vec3(p) / vec3(size - 1);
+ const auto dstVal = dstFromSrc01(src01);
+ data.at(i) = dstVal;
+ });
+ }
+
+ // -
+
+ /// p: [0, N-1] (clamps)
+ vec3 Fetch(ivec3 p) const {
+ const auto i = Index(p);
+ return data.at(i);
+ }
+
+ /// in01: [0.0, 1.0] (clamps)
+ vec3 Sample(vec3 in01) const;
+};
+
+// -
+
+/**
+Naively, it would be ideal to map directly from ycbcr to rgb,
+but headroom and footroom are problematic: For e.g. narrow-range-8-bit,
+our naive LUT would start at absolute y=0/255. However, values only start
+at y=16/255, and depending on where your first LUT sample is, you might get
+very poor approximations for y=16/255.
+Further, even for full-range-8-bit, y=-0.5 is encoded as 1/255. U and v
+aren't *as* important as y, but we should try be accurate for the min and
+max values. Additionally, it would be embarassing to get whites/greys wrong,
+so preserving u=0.0 should also be a goal.
+Finally, when using non-linear transfer functions, the linear approximation of a
+point between two samples will be fairly inaccurate.
+We preserve min and max by choosing our input range such that min and max are
+the endpoints of their LUT axis.
+We preserve accuracy (at and around) mid by choosing odd sizes for dimentions.
+
+But also, the LUT is surprisingly robust, so check if the simple version works
+before adding complexity!
+**/
+
+struct ColorspaceTransform final {
+ ColorspaceDesc srcSpace;
+ ColorspaceDesc dstSpace;
+ mat4 srcRgbTfFromSrc;
+ std::optional<PiecewiseGammaDesc> srcTf;
+ mat3 dstRgbLinFromSrcRgbLin;
+ std::optional<PiecewiseGammaDesc> dstTf;
+ mat4 dstFromDstRgbTf;
+
+ static ColorspaceTransform Create(const ColorspaceDesc& src,
+ const ColorspaceDesc& dst);
+
+ // -
+
+ vec3 DstFromSrc(vec3 src) const;
+
+ std::optional<mat4> ToMat4() const;
+
+ Lut3 ToLut3(const ivec3 size) const;
+ Lut3 ToLut3() const {
+ auto defaultSize = ivec3({31, 31, 15}); // Order of importance: G, R, B
+ if (srcSpace.yuv) {
+ defaultSize = ivec3({31, 15, 31}); // Y, Cb, Cr
+ }
+ return ToLut3(defaultSize);
+ }
+};
+
+// -
+
+struct RgbTransferTables {
+ std::vector<float> r;
+ std::vector<float> g;
+ std::vector<float> b;
+};
+float GuessGamma(const std::vector<float>& vals, float exp_guess = 1.0);
+
+static constexpr auto D65 = vec2{{0.3127, 0.3290}};
+static constexpr auto D50 = vec2{{0.34567, 0.35850}};
+mat3 XyzAFromXyzB_BradfordLinear(const vec2 xyA, const vec2 xyB);
+
+// -
+
+struct ColorProfileDesc {
+ // ICC profiles are phrased as PCS-from-encoded (PCS is CIEXYZ-D50)
+ // However, all of our colorspaces are D65, so let's normalize to that,
+ // even though it's a reversible transform.
+ color::mat4 rgbFromYcbcr = color::mat4::Identity();
+ RgbTransferTables linearFromTf;
+ color::mat3 xyzd65FromLinearRgb = color::mat3::Identity();
+
+ static ColorProfileDesc From(const ColorspaceDesc&);
+ static ColorProfileDesc From(const qcms_profile&);
+};
+
+template <class C>
+inline float SampleOutByIn(const C& outByIn, const float in) {
+ switch (outByIn.size()) {
+ case 0:
+ return in;
+ case 1:
+ return outByIn.at(0);
+ }
+ MOZ_ASSERT(outByIn.size() >= 2);
+ const auto begin = outByIn.begin();
+
+ const auto in0i = size_t(floorf(in * (outByIn.size() - 1)));
+ const auto out0_itr = begin + std::min(in0i, outByIn.size() - 2);
+
+ const auto in0 = float(out0_itr - begin) / (outByIn.size() - 1);
+ const auto out0 = *out0_itr;
+ const auto d_in = float(1) / (outByIn.size() - 1);
+ const auto d_out = *(out0_itr + 1) - *out0_itr;
+
+ const auto out = out0 + (d_out / d_in) * (in - in0);
+ // printf("SampleOutByIn(%f)->%f\n", in, out);
+ return out;
+}
+
+template <class C>
+inline float SampleInByOut(const C& outByIn, const float out) {
+ MOZ_ASSERT(outByIn.size() >= 2);
+ const auto begin = outByIn.begin();
+
+ const auto out0_itr = std::lower_bound(begin + 1, outByIn.end() - 1, out) - 1;
+
+ const auto in0 = float(out0_itr - begin) / (outByIn.size() - 1);
+ const auto out0 = *out0_itr;
+ const auto d_in = float(1) / (outByIn.size() - 1);
+ const auto d_out = *(out0_itr + 1) - *out0_itr;
+
+ // printf("%f + (%f / %f) * (%f - %f)\n", in0, d_in, d_out, out, out0);
+ const auto in = in0 + (d_in / d_out) * (out - out0);
+ // printf("SampleInByOut(%f)->%f\n", out, in);
+ return in;
+}
+
+template <class C, class FnLessEqualT = std::less_equal<typename C::value_type>>
+inline bool IsMonotonic(const C& vals, const FnLessEqualT& LessEqual = {}) {
+ bool ok = true;
+ const auto begin = vals.begin();
+ for (size_t i = 1; i < vals.size(); i++) {
+ const auto itr = begin + i;
+ ok &= LessEqual(*(itr - 1), *itr);
+ // Assert(true, [&]() {
+ // return prints("[%zu]->%f <= [%zu]->%f", i-1, *(itr-1), i, *itr);
+ // });
+ }
+ return ok;
+}
+
+template <class T, class I>
+inline std::optional<I> SeekNeq(const T& ref, const I first, const I last) {
+ const auto inc = (last - first) > 0 ? 1 : -1;
+ auto itr = first;
+ while (true) {
+ if (*itr != ref) return itr;
+ if (itr == last) return {};
+ itr += inc;
+ }
+}
+
+template <class T>
+struct TwoPoints {
+ struct {
+ T x;
+ T y;
+ } p0;
+ struct {
+ T x;
+ T y;
+ } p1;
+
+ T y(const T x) const {
+ const auto dx = p1.x - p0.x;
+ const auto dy = p1.y - p0.y;
+ return p0.y + dy / dx * (x - p0.x);
+ }
+};
+
+/// Fills `vals` with `x:[0..vals.size()-1] => line.y(x)`.
+template <class T>
+static void LinearFill(T& vals, const TwoPoints<float>& line) {
+ float x = -1;
+ for (auto& val : vals) {
+ x += 1;
+ val = line.y(x);
+ }
+}
+
+// -
+
+inline void DequantizeMonotonic(const Span<float> vals) {
+ MOZ_ASSERT(IsMonotonic(vals));
+
+ const auto first = vals.begin();
+ const auto end = vals.end();
+ if (first == end) return;
+ const auto last = end - 1;
+ if (first == last) return;
+
+ // Three monotonic cases:
+ // 1. [0,0,0,0]
+ // 2. [0,0,1,1]
+ // 3. [0,1,1,2]
+
+ const auto body_first = SeekNeq(*first, first, last);
+ if (!body_first) {
+ // E.g. [0,0,0,0]
+ return;
+ }
+
+ const auto body_last = SeekNeq(*last, last, *body_first);
+ if (!body_last) {
+ // E.g. [0,0,1,1]
+ // This isn't the most accurate, but close enough.
+ // print("#2: %s", to_str(vals).c_str());
+ LinearFill(vals, {
+ {0, *first},
+ {float(vals.size() - 1), *last},
+ });
+ // print(" -> %s\n", to_str(vals).c_str());
+ return;
+ }
+
+ // E.g. [0,1,1,2]
+ // ^^^ body
+ // => f(0.5)->0.5, f(2.5)->1.5
+ // => f(x) = f(x0) + (x-x0) * (f(x1) - f(x0)) / (x1-x0)
+ // => f(x) = f(x0) + (x-x0) * dfdx
+
+ const auto head_end = *body_first;
+ const auto head = vals.subspan(0, head_end - vals.begin());
+ const auto tail_begin = *body_last + 1;
+ const auto tail = vals.subspan(tail_begin - vals.begin());
+ // print("head tail: %s %s\n",
+ // to_str(head).c_str(),
+ // to_str(tail).c_str());
+
+ // const auto body = vals->subspan(head.size(), vals->size()-tail.size());
+ auto next_part_first = head_end;
+ while (next_part_first != tail_begin) {
+ const auto part_first = next_part_first;
+ // print("part_first: %f\n", *part_first);
+ next_part_first = *SeekNeq(*part_first, part_first, tail_begin);
+ // print("next_part_first: %f\n", *next_part_first);
+ const auto part =
+ Span<float>{part_first, size_t(next_part_first - part_first)};
+ // print("part: %s\n", to_str(part).c_str());
+ const auto prev_part_last = part_first - 1;
+ const auto part_last = next_part_first - 1;
+ const auto line = TwoPoints<float>{
+ {-0.5, (*prev_part_last + *part_first) / 2},
+ {part.size() - 0.5f, (*part_last + *next_part_first) / 2},
+ };
+ LinearFill(part, line);
+ }
+
+ static constexpr bool INFER_HEAD_TAIL_FROM_BODY_EDGE = false;
+ // Basically ignore contents of head and tail, and infer from edges of body.
+ // print("3: %s\n", to_str(vals).c_str());
+ if (!IsMonotonic(head, std::less<float>{})) {
+ if (!INFER_HEAD_TAIL_FROM_BODY_EDGE) {
+ LinearFill(head,
+ {
+ {0, *head.begin()},
+ {head.size() - 0.5f, (*(head.end() - 1) + *head_end) / 2},
+ });
+ } else {
+ LinearFill(head, {
+ {head.size() + 0.0f, *head_end},
+ {head.size() + 1.0f, *(head_end + 1)},
+ });
+ }
+ }
+ if (!IsMonotonic(tail, std::less<float>{})) {
+ if (!INFER_HEAD_TAIL_FROM_BODY_EDGE) {
+ LinearFill(tail, {
+ {-0.5, (*(tail_begin - 1) + *tail.begin()) / 2},
+ {tail.size() - 1.0f, *(tail.end() - 1)},
+ });
+ } else {
+ LinearFill(tail, {
+ {-2.0f, *(tail_begin - 2)},
+ {-1.0f, *(tail_begin - 1)},
+ });
+ }
+ }
+ // print("3: %s\n", to_str(vals).c_str());
+ MOZ_ASSERT(IsMonotonic(vals, std::less<float>{}));
+
+ // Rescale, because we tend to lose range.
+ static constexpr bool RESCALE = false;
+ if (RESCALE) {
+ const auto firstv = *first;
+ const auto lastv = *last;
+ for (auto& val : vals) {
+ val = (val - firstv) / (lastv - firstv);
+ }
+ }
+ // print("4: %s\n", to_str(vals).c_str());
+}
+
+template <class In, class Out>
+static void InvertLut(const In& lut, Out* const out_invertedLut) {
+ MOZ_ASSERT(IsMonotonic(lut));
+ auto plut = &lut;
+ auto vec = std::vector<float>{};
+ if (!IsMonotonic(lut, std::less<float>{})) {
+ // print("Not strictly monotonic...\n");
+ vec.assign(lut.begin(), lut.end());
+ DequantizeMonotonic(vec);
+ plut = &vec;
+ // print(" Now strictly monotonic: %i: %s\n",
+ // int(IsMonotonic(*plut, std::less<float>{})), to_str(*plut).c_str());
+ MOZ_ASSERT(IsMonotonic(*plut, std::less<float>{}));
+ }
+ MOZ_ASSERT(plut->size() >= 2);
+
+ auto& ret = *out_invertedLut;
+ for (size_t i_out = 0; i_out < ret.size(); i_out++) {
+ const auto f_out = i_out / float(ret.size() - 1);
+ const auto f_in = SampleInByOut(*plut, f_out);
+ ret[i_out] = f_in;
+ }
+
+ MOZ_ASSERT(IsMonotonic(ret));
+ MOZ_ASSERT(IsMonotonic(ret, std::less<float>{}));
+}
+
+// -
+
+struct ColorProfileConversionDesc {
+ // ICC profiles are phrased as PCS-from-encoded (PCS is CIEXYZ-D50)
+ color::mat4 srcRgbFromSrcYuv = color::mat4::Identity();
+ RgbTransferTables srcLinearFromSrcTf;
+ color::mat3 dstLinearFromSrcLinear = color::mat3::Identity();
+ RgbTransferTables dstTfFromDstLinear;
+
+ struct FromDesc {
+ ColorProfileDesc src;
+ ColorProfileDesc dst;
+ };
+ static ColorProfileConversionDesc From(const FromDesc&);
+
+ vec3 Apply(const vec3 src) const {
+ const auto srcRgb = vec3(srcRgbFromSrcYuv * vec4(src, 1));
+ const auto srcLinear = vec3{{
+ SampleOutByIn(srcLinearFromSrcTf.r, srcRgb.x()),
+ SampleOutByIn(srcLinearFromSrcTf.g, srcRgb.y()),
+ SampleOutByIn(srcLinearFromSrcTf.b, srcRgb.z()),
+ }};
+ const auto dstLinear = dstLinearFromSrcLinear * srcLinear;
+ const auto dstRgb = vec3{{
+ SampleOutByIn(dstTfFromDstLinear.r, dstLinear.x()),
+ SampleOutByIn(dstTfFromDstLinear.g, dstLinear.y()),
+ SampleOutByIn(dstTfFromDstLinear.b, dstLinear.z()),
+ }};
+ return dstRgb;
+ }
+};
+
+} // namespace mozilla::color
+
+#undef ASSERT
+
+#endif // MOZILLA_GFX_GL_COLORSPACES_H_