tor-browser

The Tor Browser
git clone https://git.dasho.dev/tor-browser.git
Log | Files | Refs | README | LICENSE

TestColorspaces.cpp (26755B)


      1 /* -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*- */
      2 /* This Source Code Form is subject to the terms of the Mozilla Public
      3 * License, v. 2.0. If a copy of the MPL was not distributed with this file,
      4 * You can obtain one at http://mozilla.org/MPL/2.0/. */
      5 
      6 #include "gtest/gtest.h"
      7 #include "Colorspaces.h"
      8 
      9 #include <array>
     10 #include <limits>
     11 
     12 namespace mozilla::color {
     13 mat4 YuvFromYcbcr(const YcbcrDesc&);
     14 float TfFromLinear(const PiecewiseGammaDesc&, float linear);
     15 float LinearFromTf(const PiecewiseGammaDesc&, float tf);
     16 mat3 XyzFromLinearRgb(const Chromaticities&);
     17 }  // namespace mozilla::color
     18 
     19 using namespace mozilla::color;
     20 
     21 auto Calc8From8(const ColorspaceTransform& ct, const ivec3 in8) {
     22  const auto in = vec3(in8) / vec3(255);
     23  const auto out = ct.DstFromSrc(in);
     24  const auto out8 = ivec3(round(out * vec3(255)));
     25  return out8;
     26 }
     27 
     28 auto Sample8From8(const Lut3& lut, const vec3 in8) {
     29  const auto in = in8 / vec3(255);
     30  const auto out = lut.Sample(in);
     31  const auto out8 = ivec3(round(out * vec3(255)));
     32  return out8;
     33 }
     34 
     35 TEST(Colorspaces, YcbcrDesc_Narrow8)
     36 {
     37  const auto m = YuvFromYcbcr(YcbcrDesc::Narrow8());
     38 
     39  const auto Yuv8 = [&](const ivec3 ycbcr8) {
     40    const auto ycbcr = vec4(vec3(ycbcr8) / 255, 1);
     41    const auto yuv = m * ycbcr;
     42    return ivec3(round(yuv * 255));
     43  };
     44 
     45  EXPECT_EQ(Yuv8({{16, 128, 128}}), (ivec3{{0, 0, 0}}));
     46  EXPECT_EQ(Yuv8({{17, 128, 128}}), (ivec3{{1, 0, 0}}));
     47  // y = 0.5 => (16 + 235) / 2 = 125.5
     48  EXPECT_EQ(Yuv8({{125, 128, 128}}), (ivec3{{127, 0, 0}}));
     49  EXPECT_EQ(Yuv8({{126, 128, 128}}), (ivec3{{128, 0, 0}}));
     50  EXPECT_EQ(Yuv8({{234, 128, 128}}), (ivec3{{254, 0, 0}}));
     51  EXPECT_EQ(Yuv8({{235, 128, 128}}), (ivec3{{255, 0, 0}}));
     52 
     53  // Check that we get the naive out-of-bounds behavior we'd expect:
     54  EXPECT_EQ(Yuv8({{15, 128, 128}}), (ivec3{{-1, 0, 0}}));
     55  EXPECT_EQ(Yuv8({{236, 128, 128}}), (ivec3{{256, 0, 0}}));
     56 }
     57 
     58 TEST(Colorspaces, YcbcrDesc_Full8)
     59 {
     60  const auto m = YuvFromYcbcr(YcbcrDesc::Full8());
     61 
     62  const auto Yuv8 = [&](const ivec3 ycbcr8) {
     63    const auto ycbcr = vec4(vec3(ycbcr8) / 255, 1);
     64    const auto yuv = m * ycbcr;
     65    return ivec3(round(yuv * 255));
     66  };
     67 
     68  EXPECT_EQ(Yuv8({{0, 128, 128}}), (ivec3{{0, 0, 0}}));
     69  EXPECT_EQ(Yuv8({{1, 128, 128}}), (ivec3{{1, 0, 0}}));
     70  EXPECT_EQ(Yuv8({{127, 128, 128}}), (ivec3{{127, 0, 0}}));
     71  EXPECT_EQ(Yuv8({{128, 128, 128}}), (ivec3{{128, 0, 0}}));
     72  EXPECT_EQ(Yuv8({{254, 128, 128}}), (ivec3{{254, 0, 0}}));
     73  EXPECT_EQ(Yuv8({{255, 128, 128}}), (ivec3{{255, 0, 0}}));
     74 }
     75 
     76 TEST(Colorspaces, YcbcrDesc_Float)
     77 {
     78  const auto m = YuvFromYcbcr(YcbcrDesc::Float());
     79 
     80  const auto Yuv8 = [&](const vec3 ycbcr8) {
     81    const auto ycbcr = vec4(vec3(ycbcr8) / 255, 1);
     82    const auto yuv = m * ycbcr;
     83    return ivec3(round(yuv * 255));
     84  };
     85 
     86  EXPECT_EQ(Yuv8({{0, 0.5 * 255, 0.5 * 255}}), (ivec3{{0, 0, 0}}));
     87  EXPECT_EQ(Yuv8({{1, 0.5 * 255, 0.5 * 255}}), (ivec3{{1, 0, 0}}));
     88  EXPECT_EQ(Yuv8({{127, 0.5 * 255, 0.5 * 255}}), (ivec3{{127, 0, 0}}));
     89  EXPECT_EQ(Yuv8({{128, 0.5 * 255, 0.5 * 255}}), (ivec3{{128, 0, 0}}));
     90  EXPECT_EQ(Yuv8({{254, 0.5 * 255, 0.5 * 255}}), (ivec3{{254, 0, 0}}));
     91  EXPECT_EQ(Yuv8({{255, 0.5 * 255, 0.5 * 255}}), (ivec3{{255, 0, 0}}));
     92 }
     93 
     94 TEST(Colorspaces, ColorspaceTransform_Rec709Narrow)
     95 {
     96  const auto src = ColorspaceDesc{
     97      Chromaticities::Rec709(),
     98      PiecewiseGammaDesc::Rec709(),
     99      {{YuvLumaCoeffs::Rec709(), YcbcrDesc::Narrow8()}},
    100  };
    101  const auto dst = ColorspaceDesc{
    102      Chromaticities::Rec709(),
    103      PiecewiseGammaDesc::Rec709(),
    104      {},
    105  };
    106  const auto ct = ColorspaceTransform::Create(src, dst);
    107 
    108  EXPECT_EQ(Calc8From8(ct, {{16, 128, 128}}), (ivec3{0}));
    109  EXPECT_EQ(Calc8From8(ct, {{17, 128, 128}}), (ivec3{1}));
    110  EXPECT_EQ(Calc8From8(ct, {{126, 128, 128}}), (ivec3{128}));
    111  EXPECT_EQ(Calc8From8(ct, {{234, 128, 128}}), (ivec3{254}));
    112  EXPECT_EQ(Calc8From8(ct, {{235, 128, 128}}), (ivec3{255}));
    113 
    114  // Check that we get the naive out-of-bounds behavior we'd expect:
    115  EXPECT_EQ(Calc8From8(ct, {{15, 128, 128}}), (ivec3{-1}));
    116  EXPECT_EQ(Calc8From8(ct, {{236, 128, 128}}), (ivec3{256}));
    117 }
    118 
    119 TEST(Colorspaces, LutSample_Rec709Float)
    120 {
    121  const auto src = ColorspaceDesc{
    122      Chromaticities::Rec709(),
    123      PiecewiseGammaDesc::Rec709(),
    124      {{YuvLumaCoeffs::Rec709(), YcbcrDesc::Float()}},
    125  };
    126  const auto dst = ColorspaceDesc{
    127      Chromaticities::Rec709(),
    128      PiecewiseGammaDesc::Rec709(),
    129      {},
    130  };
    131  const auto lut = ColorspaceTransform::Create(src, dst).ToLut3();
    132 
    133  EXPECT_EQ(Sample8From8(lut, {{0, 0.5 * 255, 0.5 * 255}}), (ivec3{0}));
    134  EXPECT_EQ(Sample8From8(lut, {{1, 0.5 * 255, 0.5 * 255}}), (ivec3{1}));
    135  EXPECT_EQ(Sample8From8(lut, {{127, 0.5 * 255, 0.5 * 255}}), (ivec3{127}));
    136  EXPECT_EQ(Sample8From8(lut, {{128, 0.5 * 255, 0.5 * 255}}), (ivec3{128}));
    137  EXPECT_EQ(Sample8From8(lut, {{254, 0.5 * 255, 0.5 * 255}}), (ivec3{254}));
    138  EXPECT_EQ(Sample8From8(lut, {{255, 0.5 * 255, 0.5 * 255}}), (ivec3{255}));
    139 }
    140 
    141 TEST(Colorspaces, LutSample_Rec709Narrow)
    142 {
    143  const auto src = ColorspaceDesc{
    144      Chromaticities::Rec709(),
    145      PiecewiseGammaDesc::Rec709(),
    146      {{YuvLumaCoeffs::Rec709(), YcbcrDesc::Narrow8()}},
    147  };
    148  const auto dst = ColorspaceDesc{
    149      Chromaticities::Rec709(),
    150      PiecewiseGammaDesc::Rec709(),
    151      {},
    152  };
    153  const auto lut = ColorspaceTransform::Create(src, dst).ToLut3();
    154 
    155  EXPECT_EQ(Sample8From8(lut, {{16, 128, 128}}), (ivec3{0}));
    156  EXPECT_EQ(Sample8From8(lut, {{17, 128, 128}}), (ivec3{1}));
    157  EXPECT_EQ(Sample8From8(lut, {{int((235 + 16) / 2), 128, 128}}), (ivec3{127}));
    158  EXPECT_EQ(Sample8From8(lut, {{int((235 + 16) / 2) + 1, 128, 128}}),
    159            (ivec3{128}));
    160  EXPECT_EQ(Sample8From8(lut, {{234, 128, 128}}), (ivec3{254}));
    161  EXPECT_EQ(Sample8From8(lut, {{235, 128, 128}}), (ivec3{255}));
    162 }
    163 
    164 TEST(Colorspaces, LutSample_Rec709Full)
    165 {
    166  const auto src = ColorspaceDesc{
    167      Chromaticities::Rec709(),
    168      PiecewiseGammaDesc::Rec709(),
    169      {{YuvLumaCoeffs::Rec709(), YcbcrDesc::Full8()}},
    170  };
    171  const auto dst = ColorspaceDesc{
    172      Chromaticities::Rec709(),
    173      PiecewiseGammaDesc::Rec709(),
    174      {},
    175  };
    176  const auto lut = ColorspaceTransform::Create(src, dst).ToLut3();
    177 
    178  EXPECT_EQ(Sample8From8(lut, {{0, 128, 128}}), (ivec3{0}));
    179  EXPECT_EQ(Sample8From8(lut, {{1, 128, 128}}), (ivec3{1}));
    180  EXPECT_EQ(Sample8From8(lut, {{16, 128, 128}}), (ivec3{16}));
    181  EXPECT_EQ(Sample8From8(lut, {{128, 128, 128}}), (ivec3{128}));
    182  EXPECT_EQ(Sample8From8(lut, {{235, 128, 128}}), (ivec3{235}));
    183  EXPECT_EQ(Sample8From8(lut, {{254, 128, 128}}), (ivec3{254}));
    184  EXPECT_EQ(Sample8From8(lut, {{255, 128, 128}}), (ivec3{255}));
    185 }
    186 
    187 TEST(Colorspaces, PiecewiseGammaDesc_Srgb)
    188 {
    189  const auto tf = PiecewiseGammaDesc::Srgb();
    190 
    191  EXPECT_EQ(int(roundf(TfFromLinear(tf, 0x00 / 255.0) * 255)), 0x00);
    192  EXPECT_EQ(int(roundf(TfFromLinear(tf, 0x01 / 255.0) * 255)), 0x0d);
    193  EXPECT_EQ(int(roundf(TfFromLinear(tf, 0x37 / 255.0) * 255)), 0x80);
    194  EXPECT_EQ(int(roundf(TfFromLinear(tf, 0x80 / 255.0) * 255)), 0xbc);
    195  EXPECT_EQ(int(roundf(TfFromLinear(tf, 0xfd / 255.0) * 255)), 0xfe);
    196  EXPECT_EQ(int(roundf(TfFromLinear(tf, 0xfe / 255.0) * 255)), 0xff);
    197  EXPECT_EQ(int(roundf(TfFromLinear(tf, 0xff / 255.0) * 255)), 0xff);
    198 
    199  EXPECT_EQ(int(roundf(LinearFromTf(tf, 0x00 / 255.0) * 255)), 0x00);
    200  EXPECT_EQ(int(roundf(LinearFromTf(tf, 0x01 / 255.0) * 255)), 0x00);
    201  EXPECT_EQ(int(roundf(LinearFromTf(tf, 0x06 / 255.0) * 255)), 0x00);
    202  EXPECT_EQ(int(roundf(LinearFromTf(tf, 0x07 / 255.0) * 255)), 0x01);
    203  EXPECT_EQ(int(roundf(LinearFromTf(tf, 0x0d / 255.0) * 255)), 0x01);
    204  EXPECT_EQ(int(roundf(LinearFromTf(tf, 0x80 / 255.0) * 255)), 0x37);
    205  EXPECT_EQ(int(roundf(LinearFromTf(tf, 0xbc / 255.0) * 255)), 0x80);
    206  EXPECT_EQ(int(roundf(LinearFromTf(tf, 0xfe / 255.0) * 255)), 0xfd);
    207  EXPECT_EQ(int(roundf(LinearFromTf(tf, 0xff / 255.0) * 255)), 0xff);
    208 }
    209 
    210 TEST(Colorspaces, PiecewiseGammaDesc_Rec709)
    211 {
    212  const auto tf = PiecewiseGammaDesc::Rec709();
    213 
    214  EXPECT_EQ(int(roundf(TfFromLinear(tf, 0x00 / 255.0) * 255)), 0x00);
    215  EXPECT_EQ(int(roundf(TfFromLinear(tf, 0x01 / 255.0) * 255)), 0x05);
    216  EXPECT_EQ(int(roundf(TfFromLinear(tf, 0x43 / 255.0) * 255)), 0x80);
    217  EXPECT_EQ(int(roundf(TfFromLinear(tf, 0x80 / 255.0) * 255)), 0xb4);
    218  EXPECT_EQ(int(roundf(TfFromLinear(tf, 0xfd / 255.0) * 255)), 0xfe);
    219  EXPECT_EQ(int(roundf(TfFromLinear(tf, 0xfe / 255.0) * 255)), 0xff);
    220  EXPECT_EQ(int(roundf(TfFromLinear(tf, 0xff / 255.0) * 255)), 0xff);
    221 
    222  EXPECT_EQ(int(roundf(LinearFromTf(tf, 0x00 / 255.0) * 255)), 0x00);
    223  EXPECT_EQ(int(roundf(LinearFromTf(tf, 0x01 / 255.0) * 255)), 0x00);
    224  EXPECT_EQ(int(roundf(LinearFromTf(tf, 0x02 / 255.0) * 255)), 0x00);
    225  EXPECT_EQ(int(roundf(LinearFromTf(tf, 0x03 / 255.0) * 255)), 0x01);
    226  EXPECT_EQ(int(roundf(LinearFromTf(tf, 0x05 / 255.0) * 255)), 0x01);
    227  EXPECT_EQ(int(roundf(LinearFromTf(tf, 0x80 / 255.0) * 255)), 0x43);
    228  EXPECT_EQ(int(roundf(LinearFromTf(tf, 0xb4 / 255.0) * 255)), 0x80);
    229  EXPECT_EQ(int(roundf(LinearFromTf(tf, 0xfe / 255.0) * 255)), 0xfd);
    230  EXPECT_EQ(int(roundf(LinearFromTf(tf, 0xff / 255.0) * 255)), 0xff);
    231 }
    232 
    233 TEST(Colorspaces, ColorspaceTransform_PiecewiseGammaDesc)
    234 {
    235  const auto src = ColorspaceDesc{
    236      Chromaticities::Srgb(),
    237      {},
    238      {},
    239  };
    240  const auto dst = ColorspaceDesc{
    241      Chromaticities::Srgb(),
    242      PiecewiseGammaDesc::Srgb(),
    243      {},
    244  };
    245  const auto toGamma = ColorspaceTransform::Create(src, dst);
    246  const auto toLinear = ColorspaceTransform::Create(dst, src);
    247 
    248  EXPECT_EQ(Calc8From8(toGamma, ivec3{0x00}), (ivec3{0x00}));
    249  EXPECT_EQ(Calc8From8(toGamma, ivec3{0x01}), (ivec3{0x0d}));
    250  EXPECT_EQ(Calc8From8(toGamma, ivec3{0x37}), (ivec3{0x80}));
    251  EXPECT_EQ(Calc8From8(toGamma, ivec3{0x80}), (ivec3{0xbc}));
    252  EXPECT_EQ(Calc8From8(toGamma, ivec3{0xfd}), (ivec3{0xfe}));
    253  EXPECT_EQ(Calc8From8(toGamma, ivec3{0xff}), (ivec3{0xff}));
    254 
    255  EXPECT_EQ(Calc8From8(toLinear, ivec3{0x00}), (ivec3{0x00}));
    256  EXPECT_EQ(Calc8From8(toLinear, ivec3{0x0d}), (ivec3{0x01}));
    257  EXPECT_EQ(Calc8From8(toLinear, ivec3{0x80}), (ivec3{0x37}));
    258  EXPECT_EQ(Calc8From8(toLinear, ivec3{0xbc}), (ivec3{0x80}));
    259  EXPECT_EQ(Calc8From8(toLinear, ivec3{0xfe}), (ivec3{0xfd}));
    260  EXPECT_EQ(Calc8From8(toLinear, ivec3{0xff}), (ivec3{0xff}));
    261 }
    262 
    263 // -
    264 // Actual end-to-end tests
    265 
    266 TEST(Colorspaces, SrgbFromRec709)
    267 {
    268  const auto src = ColorspaceDesc{
    269      Chromaticities::Rec709(),
    270      PiecewiseGammaDesc::Rec709(),
    271      {{YuvLumaCoeffs::Rec709(), YcbcrDesc::Narrow8()}},
    272  };
    273  const auto dst = ColorspaceDesc{
    274      Chromaticities::Srgb(),
    275      PiecewiseGammaDesc::Srgb(),
    276      {},
    277  };
    278  const auto ct = ColorspaceTransform::Create(src, dst);
    279 
    280  EXPECT_EQ(Calc8From8(ct, ivec3{{16, 128, 128}}), (ivec3{0}));
    281  EXPECT_EQ(Calc8From8(ct, ivec3{{17, 128, 128}}), (ivec3{3}));
    282  EXPECT_EQ(Calc8From8(ct, ivec3{{115, 128, 128}}), (ivec3{128}));
    283  EXPECT_EQ(Calc8From8(ct, ivec3{{126, 128, 128}}), (ivec3{140}));
    284  EXPECT_EQ(Calc8From8(ct, ivec3{{234, 128, 128}}), (ivec3{254}));
    285  EXPECT_EQ(Calc8From8(ct, ivec3{{235, 128, 128}}), (ivec3{255}));
    286 }
    287 
    288 TEST(Colorspaces, SrgbFromDisplayP3)
    289 {
    290  const auto p3C = ColorspaceDesc{
    291      Chromaticities::DisplayP3(),
    292      PiecewiseGammaDesc::DisplayP3(),
    293  };
    294  const auto srgbC = ColorspaceDesc{
    295      Chromaticities::Srgb(),
    296      PiecewiseGammaDesc::Srgb(),
    297  };
    298  const auto srgbLinearC = ColorspaceDesc{
    299      Chromaticities::Srgb(),
    300      {},
    301  };
    302  const auto srgbFromP3 = ColorspaceTransform::Create(p3C, srgbC);
    303  const auto srgbLinearFromP3 = ColorspaceTransform::Create(p3C, srgbLinearC);
    304 
    305  // E.g.
    306  // https://colorjs.io/apps/convert/?color=color(display-p3%200.4%200.8%200.4)&precision=4
    307  auto srgb = srgbFromP3.DstFromSrc(vec3{{0.4, 0.8, 0.4}});
    308  EXPECT_NEAR(srgb.x(), 0.179, 0.001);
    309  EXPECT_NEAR(srgb.y(), 0.812, 0.001);
    310  EXPECT_NEAR(srgb.z(), 0.342, 0.001);
    311  auto srgbLinear = srgbLinearFromP3.DstFromSrc(vec3{{0.4, 0.8, 0.4}});
    312  EXPECT_NEAR(srgbLinear.x(), 0.027, 0.001);
    313  EXPECT_NEAR(srgbLinear.y(), 0.624, 0.001);
    314  EXPECT_NEAR(srgbLinear.z(), 0.096, 0.001);
    315 }
    316 
    317 TEST(Colorspaces, DisplayP3FromSrgb)
    318 {
    319  const auto p3C = ColorspaceDesc{
    320      Chromaticities::DisplayP3(),
    321      PiecewiseGammaDesc::DisplayP3(),
    322  };
    323  const auto srgbC = ColorspaceDesc{
    324      Chromaticities::Srgb(),
    325      PiecewiseGammaDesc::Srgb(),
    326  };
    327  const auto p3FromSrgb = ColorspaceTransform::Create(srgbC, p3C);
    328 
    329  // E.g.
    330  // https://colorjs.io/apps/convert/?color=color(srgb%200.4%200.8%200.4)&precision=4
    331  auto p3 = p3FromSrgb.DstFromSrc(vec3{{0.4, 0.8, 0.4}});
    332  EXPECT_NEAR(p3.x(), 0.502, 0.001);
    333  EXPECT_NEAR(p3.y(), 0.791, 0.001);
    334  EXPECT_NEAR(p3.z(), 0.445, 0.001);
    335 }
    336 
    337 // -
    338 
    339 template <class Fn, class Tuple, size_t... I>
    340 constexpr auto map_tups_seq(const Tuple& a, const Tuple& b, const Fn& fn,
    341                            std::index_sequence<I...>) {
    342  return std::tuple{fn(std::get<I>(a), std::get<I>(b))...};
    343 }
    344 template <class Fn, class Tuple>
    345 constexpr auto map_tups(const Tuple& a, const Tuple& b, const Fn& fn) {
    346  return map_tups_seq(a, b, fn,
    347                      std::make_index_sequence<std::tuple_size_v<Tuple>>{});
    348 }
    349 
    350 template <class Fn, class Tuple>
    351 constexpr auto cmp_tups_all(const Tuple& a, const Tuple& b, const Fn& fn) {
    352  bool all = true;
    353  map_tups(a, b, [&](const auto& a, const auto& b) { return all &= fn(a, b); });
    354  return all;
    355 }
    356 
    357 struct Stats {
    358  double mean = 0;
    359  double variance = 0;
    360  double min = std::numeric_limits<double>::infinity();
    361  double max = -std::numeric_limits<double>::infinity();
    362 
    363  template <class T>
    364  static Stats For(const T& iterable) {
    365    auto ret = Stats{};
    366    for (const auto& cur : iterable) {
    367      ret.mean += cur;
    368      ret.min = std::min(ret.min, cur);
    369      ret.max = std::max(ret.max, cur);
    370    }
    371    ret.mean /= iterable.size();
    372    // Gather mean first before we can calc variance.
    373    for (const auto& cur : iterable) {
    374      ret.variance += pow(cur - ret.mean, 2);
    375    }
    376    ret.variance /= iterable.size();
    377    return ret;
    378  }
    379 
    380  template <class T, class U>
    381  static Stats Diff(const T& a, const U& b) {
    382    MOZ_ASSERT(a.size() == b.size());
    383    std::vector<double> diff;
    384    diff.reserve(a.size());
    385    for (size_t i = 0; i < diff.capacity(); i++) {
    386      diff.push_back(a[i] - b[i]);
    387    }
    388    return Stats::For(diff);
    389  }
    390 
    391  double standardDeviation() const { return sqrt(variance); }
    392 
    393  friend std::ostream& operator<<(std::ostream& s, const Stats& a) {
    394    return s << "Stats"
    395             << "{ mean:" << a.mean << ", stddev:" << a.standardDeviation()
    396             << ", min:" << a.min << ", max:" << a.max << " }";
    397  }
    398 
    399  struct Error {
    400    double absmean = std::numeric_limits<double>::infinity();
    401    double stddev = std::numeric_limits<double>::infinity();
    402    double absmax = std::numeric_limits<double>::infinity();
    403 
    404    constexpr auto Fields() const { return std::tie(absmean, stddev, absmax); }
    405 
    406    template <class Fn>
    407    friend constexpr bool cmp_all(const Error& a, const Error& b,
    408                                  const Fn& fn) {
    409      return cmp_tups_all(a.Fields(), b.Fields(), fn);
    410    }
    411    friend constexpr bool operator<(const Error& a, const Error& b) {
    412      return cmp_all(a, b, [](const auto& a, const auto& b) { return a < b; });
    413    }
    414    friend constexpr bool operator<=(const Error& a, const Error& b) {
    415      return cmp_all(a, b, [](const auto& a, const auto& b) { return a <= b; });
    416    }
    417 
    418    friend std::ostream& operator<<(std::ostream& s, const Error& a) {
    419      return s << "Stats::Error"
    420               << "{ absmean:" << a.absmean << ", stddev:" << a.stddev
    421               << ", absmax:" << a.absmax << " }";
    422    }
    423  };
    424 
    425  operator Error() const {
    426    return {abs(mean), standardDeviation(), std::max(abs(min), abs(max))};
    427  }
    428 };
    429 static_assert(Stats::Error{0, 0, 0} < Stats::Error{1, 1, 1});
    430 static_assert(!(Stats::Error{0, 1, 0} < Stats::Error{1, 1, 1}));
    431 static_assert(Stats::Error{0, 1, 0} <= Stats::Error{1, 1, 1});
    432 static_assert(!(Stats::Error{0, 2, 0} <= Stats::Error{1, 1, 1}));
    433 
    434 // -
    435 
    436 static Stats StatsForLutError(const ColorspaceTransform& ct,
    437                              const ivec3 srcQuants, const ivec3 dstQuants) {
    438  const auto lut = ct.ToLut3();
    439 
    440  const auto dstScale = vec3(dstQuants - 1);
    441 
    442  std::vector<double> quantErrors;
    443  quantErrors.reserve(srcQuants.x() * srcQuants.y() * srcQuants.z());
    444  ForEachSampleWithin(srcQuants, [&](const vec3& src) {
    445    const auto sampled = lut.Sample(src);
    446    const auto actual = ct.DstFromSrc(src);
    447    const auto isampled = ivec3(round(sampled * dstScale));
    448    const auto iactual = ivec3(round(actual * dstScale));
    449    const auto ierr = abs(isampled - iactual);
    450    const auto quantError = dot(ierr, ivec3{1});
    451    quantErrors.push_back(quantError);
    452    if (quantErrors.size() % 100000 == 0) {
    453      printf("%zu of %zu\n", quantErrors.size(), quantErrors.capacity());
    454    }
    455  });
    456 
    457  const auto quantErrStats = Stats::For(quantErrors);
    458  return quantErrStats;
    459 }
    460 
    461 TEST(Colorspaces, LutError_Rec709Full_Rec709Rgb)
    462 {
    463  const auto src = ColorspaceDesc{
    464      Chromaticities::Rec709(),
    465      PiecewiseGammaDesc::Rec709(),
    466      {{YuvLumaCoeffs::Rec709(), YcbcrDesc::Full8()}},
    467  };
    468  const auto dst = ColorspaceDesc{
    469      Chromaticities::Rec709(),
    470      PiecewiseGammaDesc::Rec709(),
    471      {},
    472  };
    473  const auto ct = ColorspaceTransform::Create(src, dst);
    474  const auto stats = StatsForLutError(ct, ivec3{64}, ivec3{256});
    475  EXPECT_NEAR(stats.mean, 0.000, 0.001);
    476  EXPECT_NEAR(stats.standardDeviation(), 0.008, 0.001);
    477  EXPECT_NEAR(stats.min, 0, 0.001);
    478  EXPECT_NEAR(stats.max, 1, 0.001);
    479 }
    480 
    481 TEST(Colorspaces, LutError_Rec709Full_Srgb)
    482 {
    483  const auto src = ColorspaceDesc{
    484      Chromaticities::Rec709(),
    485      PiecewiseGammaDesc::Rec709(),
    486      {{YuvLumaCoeffs::Rec709(), YcbcrDesc::Full8()}},
    487  };
    488  const auto dst = ColorspaceDesc{
    489      Chromaticities::Srgb(),
    490      PiecewiseGammaDesc::Srgb(),
    491      {},
    492  };
    493  const auto ct = ColorspaceTransform::Create(src, dst);
    494  const auto stats = StatsForLutError(ct, ivec3{64}, ivec3{256});
    495  EXPECT_NEAR(stats.mean, 0.530, 0.001);
    496  EXPECT_NEAR(stats.standardDeviation(), 1.674, 0.001);
    497  EXPECT_NEAR(stats.min, 0, 0.001);
    498  EXPECT_NEAR(stats.max, 17, 0.001);
    499 }
    500 
    501 // -
    502 // https://www.reedbeta.com/blog/python-like-enumerate-in-cpp17/
    503 
    504 template <typename T, typename TIter = decltype(std::begin(std::declval<T>())),
    505          typename = decltype(std::end(std::declval<T>()))>
    506 constexpr auto enumerate(T&& iterable) {
    507  struct iterator {
    508    size_t i;
    509    TIter iter;
    510    bool operator!=(const iterator& other) const { return iter != other.iter; }
    511    void operator++() {
    512      ++i;
    513      ++iter;
    514    }
    515    auto operator*() const { return std::tie(i, *iter); }
    516  };
    517  struct iterable_wrapper {
    518    T iterable;
    519    auto begin() { return iterator{0, std::begin(iterable)}; }
    520    auto end() { return iterator{0, std::end(iterable)}; }
    521  };
    522  return iterable_wrapper{std::forward<T>(iterable)};
    523 }
    524 
    525 inline auto MakeLinear(const float from, const float to, const int n) {
    526  std::vector<float> ret;
    527  ret.resize(n);
    528  for (auto [i, val] : enumerate(ret)) {
    529    const auto t = i / float(ret.size() - 1);
    530    val = from + (to - from) * t;
    531  }
    532  return ret;
    533 };
    534 
    535 inline auto MakeGamma(const float exp, const int n) {
    536  std::vector<float> ret;
    537  ret.resize(n);
    538  for (auto [i, val] : enumerate(ret)) {
    539    const auto t = i / float(ret.size() - 1);
    540    val = powf(t, exp);
    541  }
    542  return ret;
    543 };
    544 
    545 // -
    546 
    547 TEST(Colorspaces, GuessGamma)
    548 {
    549  EXPECT_NEAR(GuessGamma(MakeGamma(1, 11)), 1.0, 0);
    550  EXPECT_NEAR(GuessGamma(MakeGamma(2.2, 11)), 2.2, 4.8e-8);
    551  EXPECT_NEAR(GuessGamma(MakeGamma(1 / 2.2, 11)), 1 / 2.2, 1.7e-7);
    552 }
    553 
    554 // -
    555 
    556 template <class T, class U>
    557 float StdDev(const T& test, const U& ref) {
    558  float sum = 0;
    559  for (size_t i = 0; i < test.size(); i++) {
    560    const auto diff = test[i] - ref[i];
    561    sum += diff * diff;
    562  }
    563  const auto variance = sum / test.size();
    564  return sqrt(variance);
    565 }
    566 
    567 template <class T>
    568 inline void AutoLinearFill(T& vals) {
    569  LinearFill(vals, {
    570                       {0, 0},
    571                       {vals.size() - 1.0f, 1},
    572                   });
    573 }
    574 
    575 template <class T, class... More>
    576 auto MakeArray(const T& a0, const More&... args) {
    577  return std::array<T, 1 + sizeof...(More)>{a0, static_cast<float>(args)...};
    578 }
    579 
    580 TEST(Colorspaces, LinearFill)
    581 {
    582  EXPECT_NEAR(StdDev(MakeLinear(0, 1, 3), MakeArray<float>(0, 0.5, 1)), 0,
    583              0.001);
    584 
    585  auto vals = std::vector<float>(3);
    586  LinearFill(vals, {
    587                       {0, 0},
    588                       {vals.size() - 1.0f, 1},
    589                   });
    590  EXPECT_NEAR(StdDev(vals, MakeArray<float>(0, 0.5, 1)), 0, 0.001);
    591 
    592  LinearFill(vals, {
    593                       {0, 1},
    594                       {vals.size() - 1.0f, 0},
    595                   });
    596  EXPECT_NEAR(StdDev(vals, MakeArray<float>(1, 0.5, 0)), 0, 0.001);
    597 }
    598 
    599 TEST(Colorspaces, DequantizeMonotonic)
    600 {
    601  auto orig = std::vector<float>{0, 0, 0, 1, 1, 2};
    602  auto vals = orig;
    603  EXPECT_TRUE(IsMonotonic(vals));
    604  EXPECT_TRUE(!IsMonotonic(vals, std::less<float>{}));
    605  DequantizeMonotonic(vals);
    606  EXPECT_TRUE(IsMonotonic(vals, std::less<float>{}));
    607  EXPECT_LT(StdDev(vals, orig),
    608            StdDev(MakeLinear(orig.front(), orig.back(), vals.size()), orig));
    609 }
    610 
    611 TEST(Colorspaces, InvertLut)
    612 {
    613  const auto linear = MakeLinear(0, 1, 256);
    614  auto linearFromSrgb = linear;
    615  for (auto& val : linearFromSrgb) {
    616    val = powf(val, 2.2);
    617  }
    618  auto srgbFromLinearExpected = linear;
    619  for (auto& val : srgbFromLinearExpected) {
    620    val = powf(val, 1 / 2.2);
    621  }
    622 
    623  auto srgbFromLinearViaInvert = linearFromSrgb;
    624  InvertLut(linearFromSrgb, &srgbFromLinearViaInvert);
    625  // I just want to appreciate that InvertLut is a non-analytical approximation,
    626  // and yet it's extraordinarily close to the analytical inverse.
    627  EXPECT_LE(Stats::Diff(srgbFromLinearViaInvert, srgbFromLinearExpected),
    628            (Stats::Error{3e-6, 3e-6, 3e-5}));
    629 
    630  const auto srcSrgb = MakeLinear(0, 1, 256);
    631  auto roundtripSrgb = srcSrgb;
    632  for (auto& srgb : roundtripSrgb) {
    633    const auto linear = SampleOutByIn(linearFromSrgb, srgb);
    634    const auto srgb2 = SampleOutByIn(srgbFromLinearViaInvert, linear);
    635    // printf("[%f] %f -> %f -> %f\n", srgb2-srgb, srgb, linear, srgb2);
    636    srgb = srgb2;
    637  }
    638  EXPECT_LE(Stats::Diff(roundtripSrgb, srcSrgb),
    639            (Stats::Error{0.0013, 0.0046, 0.023}));
    640 }
    641 
    642 TEST(Colorspaces, XyzFromLinearRgb)
    643 {
    644  const auto xyzd65FromLinearRgb = XyzFromLinearRgb(Chromaticities::Srgb());
    645 
    646  // http://www.brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html
    647  const auto XYZD65_FROM_LINEAR_RGB = mat3({
    648      vec3{{0.4124564, 0.3575761, 0.1804375}},
    649      vec3{{0.2126729, 0.7151522, 0.0721750}},
    650      vec3{{0.0193339, 0.1191920, 0.9503041}},
    651  });
    652  EXPECT_NEAR(sqrt(dotDifference(xyzd65FromLinearRgb, XYZD65_FROM_LINEAR_RGB)),
    653              0, 0.001);
    654 }
    655 
    656 TEST(Colorspaces, ColorProfileConversionDesc_SrgbFromRec709)
    657 {
    658  const auto srgb = ColorProfileDesc::From({
    659      Chromaticities::Srgb(),
    660      PiecewiseGammaDesc::Srgb(),
    661  });
    662  const auto rec709 = ColorProfileDesc::From({
    663      Chromaticities::Rec709(),
    664      PiecewiseGammaDesc::Rec709(),
    665  });
    666 
    667  {
    668    const auto conv = ColorProfileConversionDesc::From({
    669        .src = srgb,
    670        .dst = srgb,
    671    });
    672    auto src = vec3(16.0);
    673    auto dst = conv.DstFromSrc(src / 255) * 255;
    674 
    675    const auto tfa = PiecewiseGammaDesc::Srgb();
    676    const auto tfb = PiecewiseGammaDesc::Srgb();
    677    const auto expected =
    678        TfFromLinear(tfb, LinearFromTf(tfa, src.x() / 255)) * 255;
    679 
    680    printf("%f %f %f\n", src.x(), src.y(), src.z());
    681    printf("%f %f %f\n", dst.x(), dst.y(), dst.z());
    682    EXPECT_LT(Stats::Diff(dst.data, vec3(expected).data), (Stats::Error{0.42}));
    683  }
    684  {
    685    const auto conv = ColorProfileConversionDesc::From({
    686        .src = rec709,
    687        .dst = rec709,
    688    });
    689    auto src = vec3(16.0);
    690    auto dst = conv.DstFromSrc(src / 255) * 255;
    691 
    692    const auto tfa = PiecewiseGammaDesc::Rec709();
    693    const auto tfb = PiecewiseGammaDesc::Rec709();
    694    const auto expected =
    695        TfFromLinear(tfb, LinearFromTf(tfa, src.x() / 255)) * 255;
    696 
    697    printf("%f %f %f\n", src.x(), src.y(), src.z());
    698    printf("%f %f %f\n", dst.x(), dst.y(), dst.z());
    699    EXPECT_LT(Stats::Diff(dst.data, vec3(expected).data), (Stats::Error{1e-6}));
    700  }
    701  {
    702    const auto conv = ColorProfileConversionDesc::From({
    703        .src = rec709,
    704        .dst = srgb,
    705    });
    706    auto src = vec3(16.0);
    707    auto dst = conv.DstFromSrc(src / 255) * 255;
    708 
    709    const auto tfa = PiecewiseGammaDesc::Rec709();
    710    const auto tfb = PiecewiseGammaDesc::Srgb();
    711    const auto expected =
    712        TfFromLinear(tfb, LinearFromTf(tfa, src.x() / 255)) * 255;
    713    printf("expected: %f\n", expected);
    714    printf("%f %f %f\n", src.x(), src.y(), src.z());
    715    printf("%f %f %f\n", dst.x(), dst.y(), dst.z());
    716    EXPECT_LT(Stats::Diff(dst.data, vec3(expected).data), (Stats::Error{0.12}));
    717  }
    718 }
    719 
    720 TEST(Colorspaces, SampleOutByIn_NegativeInputs)
    721 {
    722  const auto tf = MakeGamma(1.0 / 2.2, 256);
    723  EXPECT_LT(SampleOutByIn(tf, -0.1f), 0.0f);
    724 }
    725 
    726 TEST(Colorspaces, ColorProfileConversionDesc_Srgb_Displayp3)
    727 {
    728  const auto srgb = ColorProfileDesc::From({
    729      Chromaticities::Srgb(),
    730      PiecewiseGammaDesc::Srgb(),
    731  });
    732  const auto displayp3 = ColorProfileDesc::From({
    733      Chromaticities::DisplayP3(),
    734      PiecewiseGammaDesc::DisplayP3(),
    735  });
    736  const auto srgbToDisplayp3 = ColorProfileConversionDesc::From({
    737      .src = srgb,
    738      .dst = displayp3,
    739  });
    740  const auto displayp3ToSrgb = ColorProfileConversionDesc::From({
    741      .src = displayp3,
    742      .dst = srgb,
    743  });
    744 
    745  const auto VecNear = [](const vec3& was, const vec3& expected,
    746                          const vec3& max_err) {
    747    const auto err = abs(was - expected);
    748    const auto err_above_max_err = err - max_err;
    749    bool ok = true;
    750    for (const auto v : err_above_max_err.data) {
    751      ok &= v <= 0.0;
    752    }
    753    return ok;
    754  };
    755 
    756  // https://apps.colorjs.io/convert/?color=color(srgb%201.0%200%200)&precision=4
    757  EXPECT_PRED3(VecNear, srgbToDisplayp3.DstFromSrc(vec3({1.0, 0.0, 0.0})),
    758               vec3({0.9175, 0.2003, 0.1386}), vec3(0.0004));
    759  EXPECT_PRED3(VecNear, srgbToDisplayp3.DstFromSrc(vec3({0.0, 1.0, 0.0})),
    760               vec3({0.4584, 0.9853, 0.2983}), vec3(0.0001));
    761  EXPECT_PRED3(VecNear, srgbToDisplayp3.DstFromSrc(vec3({0.0, 0.0, 1.0})),
    762               vec3({0.0000, 0.0000, 0.9596}), vec3(0.0001));
    763  EXPECT_PRED3(VecNear,
    764               displayp3ToSrgb.DstFromSrc(vec3({0.9175, 0.2003, 0.1386})),
    765               vec3({1.0, 0.0, 0.0}), vec3(0.0002));
    766  EXPECT_PRED3(VecNear,
    767               displayp3ToSrgb.DstFromSrc(vec3({0.4584, 0.9853, 0.2983})),
    768               vec3({0.0, 1.0, 0.0}), vec3(0.0003));
    769  EXPECT_PRED3(VecNear,
    770               displayp3ToSrgb.DstFromSrc(vec3({0.0000, 0.0000, 0.9596})),
    771               vec3({0.0, 0.0, 1.0}), vec3(0.0001));
    772 }