Adding upstream version 115.7.0esr.upstream/115.7.0esr

Signed-off-by: Daniel Baumann <daniel.baumann@progress-linux.org>

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_