gfx/gl/gtest/TestColorspaces.cpp

   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 // -
 318
 319 template <class Fn, class Tuple, size_t... I>
 320 constexpr auto map_tups_seq(const Tuple& a, const Tuple& b, const Fn& fn,
 321                             std::index_sequence<I...>) {
 322   return std::tuple{fn(std::get<I>(a), std::get<I>(b))...};
 323 }
 324 template <class Fn, class Tuple>
 325 constexpr auto map_tups(const Tuple& a, const Tuple& b, const Fn& fn) {
 326   return map_tups_seq(a, b, fn,
 327                       std::make_index_sequence<std::tuple_size_v<Tuple>>{});
 328 }
 329
 330 template <class Fn, class Tuple>
 331 constexpr auto cmp_tups_all(const Tuple& a, const Tuple& b, const Fn& fn) {
 332   bool all = true;
 333   map_tups(a, b, [&](const auto& a, const auto& b) { return all &= fn(a, b); });
 334   return all;
 335 }
 336
 337 struct Stats {
 338   double mean = 0;
 339   double variance = 0;
 340   double min = std::numeric_limits<double>::infinity();
 341   double max = -std::numeric_limits<double>::infinity();
 342
 343   template <class T>
 344   static Stats For(const T& iterable) {
 345     auto ret = Stats{};
 346     for (const auto& cur : iterable) {
 347       ret.mean += cur;
 348       ret.min = std::min(ret.min, cur);
 349       ret.max = std::max(ret.max, cur);
 350     }
 351     ret.mean /= iterable.size();
 352     // Gather mean first before we can calc variance.
 353     for (const auto& cur : iterable) {
 354       ret.variance += pow(cur - ret.mean, 2);
 355     }
 356     ret.variance /= iterable.size();
 357     return ret;
 358   }
 359
 360   template <class T, class U>
 361   static Stats Diff(const T& a, const U& b) {
 362     MOZ_ASSERT(a.size() == b.size());
 363     std::vector<double> diff;
 364     diff.reserve(a.size());
 365     for (size_t i = 0; i < diff.capacity(); i++) {
 366       diff.push_back(a[i] - b[i]);
 367     }
 368     return Stats::For(diff);
 369   }
 370
 371   double standardDeviation() const { return sqrt(variance); }
 372
 373   friend std::ostream& operator<<(std::ostream& s, const Stats& a) {
 374     return s << "Stats"
 375              << "{ mean:" << a.mean << ", stddev:" << a.standardDeviation()
 376              << ", min:" << a.min << ", max:" << a.max << " }";
 377   }
 378
 379   struct Error {
 380     double absmean = std::numeric_limits<double>::infinity();
 381     double stddev = std::numeric_limits<double>::infinity();
 382     double absmax = std::numeric_limits<double>::infinity();
 383
 384     constexpr auto Fields() const { return std::tie(absmean, stddev, absmax); }
 385
 386     template <class Fn>
 387     friend constexpr bool cmp_all(const Error& a, const Error& b,
 388                                   const Fn& fn) {
 389       return cmp_tups_all(a.Fields(), b.Fields(), fn);
 390     }
 391     friend constexpr bool operator<(const Error& a, const Error& b) {
 392       return cmp_all(a, b, [](const auto& a, const auto& b) { return a < b; });
 393     }
 394     friend constexpr bool operator<=(const Error& a, const Error& b) {
 395       return cmp_all(a, b, [](const auto& a, const auto& b) { return a <= b; });
 396     }
 397
 398     friend std::ostream& operator<<(std::ostream& s, const Error& a) {
 399       return s << "Stats::Error"
 400                << "{ absmean:" << a.absmean << ", stddev:" << a.stddev
 401                << ", absmax:" << a.absmax << " }";
 402     }
 403   };
 404
 405   operator Error() const {
 406     return {abs(mean), standardDeviation(), std::max(abs(min), abs(max))};
 407   }
 408 };
 409 static_assert(Stats::Error{0, 0, 0} < Stats::Error{1, 1, 1});
 410 static_assert(!(Stats::Error{0, 1, 0} < Stats::Error{1, 1, 1}));
 411 static_assert(Stats::Error{0, 1, 0} <= Stats::Error{1, 1, 1});
 412 static_assert(!(Stats::Error{0, 2, 0} <= Stats::Error{1, 1, 1}));
 413
 414 // -
 415
 416 static Stats StatsForLutError(const ColorspaceTransform& ct,
 417                               const ivec3 srcQuants, const ivec3 dstQuants) {
 418   const auto lut = ct.ToLut3();
 419
 420   const auto dstScale = vec3(dstQuants - 1);
 421
 422   std::vector<double> quantErrors;
 423   quantErrors.reserve(srcQuants.x() * srcQuants.y() * srcQuants.z());
 424   ForEachSampleWithin(srcQuants, [&](const vec3& src) {
 425     const auto sampled = lut.Sample(src);
 426     const auto actual = ct.DstFromSrc(src);
 427     const auto isampled = ivec3(round(sampled * dstScale));
 428     const auto iactual = ivec3(round(actual * dstScale));
 429     const auto ierr = abs(isampled - iactual);
 430     const auto quantError = dot(ierr, ivec3{1});
 431     quantErrors.push_back(quantError);
 432     if (quantErrors.size() % 100000 == 0) {
 433       printf("%zu of %zu\n", quantErrors.size(), quantErrors.capacity());
 434     }
 435   });
 436
 437   const auto quantErrStats = Stats::For(quantErrors);
 438   return quantErrStats;
 439 }
 440
 441 TEST(Colorspaces, LutError_Rec709Full_Rec709Rgb)
 442 {
 443   const auto src = ColorspaceDesc{
 444       Chromaticities::Rec709(),
 445       PiecewiseGammaDesc::Rec709(),
 446       {{YuvLumaCoeffs::Rec709(), YcbcrDesc::Full8()}},
 447   };
 448   const auto dst = ColorspaceDesc{
 449       Chromaticities::Rec709(),
 450       PiecewiseGammaDesc::Rec709(),
 451       {},
 452   };
 453   const auto ct = ColorspaceTransform::Create(src, dst);
 454   const auto stats = StatsForLutError(ct, ivec3{64}, ivec3{256});
 455   EXPECT_NEAR(stats.mean, 0.000, 0.001);
 456   EXPECT_NEAR(stats.standardDeviation(), 0.008, 0.001);
 457   EXPECT_NEAR(stats.min, 0, 0.001);
 458   EXPECT_NEAR(stats.max, 1, 0.001);
 459 }
 460
 461 TEST(Colorspaces, LutError_Rec709Full_Srgb)
 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::Srgb(),
 470       PiecewiseGammaDesc::Srgb(),
 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.530, 0.001);
 476   EXPECT_NEAR(stats.standardDeviation(), 1.674, 0.001);
 477   EXPECT_NEAR(stats.min, 0, 0.001);
 478   EXPECT_NEAR(stats.max, 17, 0.001);
 479 }
 480
 481 // -
 482 // https://www.reedbeta.com/blog/python-like-enumerate-in-cpp17/
 483
 484 template <typename T, typename TIter = decltype(std::begin(std::declval<T>())),
 485           typename = decltype(std::end(std::declval<T>()))>
 486 constexpr auto enumerate(T&& iterable) {
 487   struct iterator {
 488     size_t i;
 489     TIter iter;
 490     bool operator!=(const iterator& other) const { return iter != other.iter; }
 491     void operator++() {
 492       ++i;
 493       ++iter;
 494     }
 495     auto operator*() const { return std::tie(i, *iter); }
 496   };
 497   struct iterable_wrapper {
 498     T iterable;
 499     auto begin() { return iterator{0, std::begin(iterable)}; }
 500     auto end() { return iterator{0, std::end(iterable)}; }
 501   };
 502   return iterable_wrapper{std::forward<T>(iterable)};
 503 }
 504
 505 inline auto MakeLinear(const float from, const float to, const int n) {
 506   std::vector<float> ret;
 507   ret.resize(n);
 508   for (auto [i, val] : enumerate(ret)) {
 509     const auto t = i / float(ret.size() - 1);
 510     val = from + (to - from) * t;
 511   }
 512   return ret;
 513 };
 514
 515 inline auto MakeGamma(const float exp, const int n) {
 516   std::vector<float> ret;
 517   ret.resize(n);
 518   for (auto [i, val] : enumerate(ret)) {
 519     const auto t = i / float(ret.size() - 1);
 520     val = powf(t, exp);
 521   }
 522   return ret;
 523 };
 524
 525 // -
 526
 527 TEST(Colorspaces, GuessGamma)
 528 {
 529   EXPECT_NEAR(GuessGamma(MakeGamma(1, 11)), 1.0, 0);
 530   EXPECT_NEAR(GuessGamma(MakeGamma(2.2, 11)), 2.2, 4.8e-8);
 531   EXPECT_NEAR(GuessGamma(MakeGamma(1 / 2.2, 11)), 1 / 2.2, 1.7e-7);
 532 }
 533
 534 // -
 535
 536 template <class T, class U>
 537 float StdDev(const T& test, const U& ref) {
 538   float sum = 0;
 539   for (size_t i = 0; i < test.size(); i++) {
 540     const auto diff = test[i] - ref[i];
 541     sum += diff * diff;
 542   }
 543   const auto variance = sum / test.size();
 544   return sqrt(variance);
 545 }
 546
 547 template <class T>
 548 inline void AutoLinearFill(T& vals) {
 549   LinearFill(vals, {
 550                        {0, 0},
 551                        {vals.size() - 1.0f, 1},
 552                    });
 553 }
 554
 555 template <class T, class... More>
 556 auto MakeArray(const T& a0, const More&... args) {
 557   return std::array<T, 1 + sizeof...(More)>{a0, static_cast<float>(args)...};
 558 }
 559
 560 TEST(Colorspaces, LinearFill)
 561 {
 562   EXPECT_NEAR(StdDev(MakeLinear(0, 1, 3), MakeArray<float>(0, 0.5, 1)), 0,
 563               0.001);
 564
 565   auto vals = std::vector<float>(3);
 566   LinearFill(vals, {
 567                        {0, 0},
 568                        {vals.size() - 1.0f, 1},
 569                    });
 570   EXPECT_NEAR(StdDev(vals, MakeArray<float>(0, 0.5, 1)), 0, 0.001);
 571
 572   LinearFill(vals, {
 573                        {0, 1},
 574                        {vals.size() - 1.0f, 0},
 575                    });
 576   EXPECT_NEAR(StdDev(vals, MakeArray<float>(1, 0.5, 0)), 0, 0.001);
 577 }
 578
 579 TEST(Colorspaces, DequantizeMonotonic)
 580 {
 581   auto orig = std::vector<float>{0, 0, 0, 1, 1, 2};
 582   auto vals = orig;
 583   EXPECT_TRUE(IsMonotonic(vals));
 584   EXPECT_TRUE(!IsMonotonic(vals, std::less<float>{}));
 585   DequantizeMonotonic(vals);
 586   EXPECT_TRUE(IsMonotonic(vals, std::less<float>{}));
 587   EXPECT_LT(StdDev(vals, orig),
 588             StdDev(MakeLinear(orig.front(), orig.back(), vals.size()), orig));
 589 }
 590
 591 TEST(Colorspaces, InvertLut)
 592 {
 593   const auto linear = MakeLinear(0, 1, 256);
 594   auto linearFromSrgb = linear;
 595   for (auto& val : linearFromSrgb) {
 596     val = powf(val, 2.2);
 597   }
 598   auto srgbFromLinearExpected = linear;
 599   for (auto& val : srgbFromLinearExpected) {
 600     val = powf(val, 1 / 2.2);
 601   }
 602
 603   auto srgbFromLinearViaInvert = linearFromSrgb;
 604   InvertLut(linearFromSrgb, &srgbFromLinearViaInvert);
 605   // I just want to appreciate that InvertLut is a non-analytical approximation,
 606   // and yet it's extraordinarily close to the analytical inverse.
 607   EXPECT_LE(Stats::Diff(srgbFromLinearViaInvert, srgbFromLinearExpected),
 608             (Stats::Error{3e-6, 3e-6, 3e-5}));
 609
 610   const auto srcSrgb = MakeLinear(0, 1, 256);
 611   auto roundtripSrgb = srcSrgb;
 612   for (auto& srgb : roundtripSrgb) {
 613     const auto linear = SampleOutByIn(linearFromSrgb, srgb);
 614     const auto srgb2 = SampleOutByIn(srgbFromLinearViaInvert, linear);
 615     // printf("[%f] %f -> %f -> %f\n", srgb2-srgb, srgb, linear, srgb2);
 616     srgb = srgb2;
 617   }
 618   EXPECT_LE(Stats::Diff(roundtripSrgb, srcSrgb),
 619             (Stats::Error{0.0013, 0.0046, 0.023}));
 620 }
 621
 622 TEST(Colorspaces, XyzFromLinearRgb)
 623 {
 624   const auto xyzd65FromLinearRgb = XyzFromLinearRgb(Chromaticities::Srgb());
 625
 626   // http://www.brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html
 627   const auto XYZD65_FROM_LINEAR_RGB = mat3({
 628       vec3{{0.4124564, 0.3575761, 0.1804375}},
 629       vec3{{0.2126729, 0.7151522, 0.0721750}},
 630       vec3{{0.0193339, 0.1191920, 0.9503041}},
 631   });
 632   EXPECT_NEAR(sqrt(dotDifference(xyzd65FromLinearRgb, XYZD65_FROM_LINEAR_RGB)),
 633               0, 0.001);
 634 }
 635
 636 TEST(Colorspaces, ColorProfileConversionDesc_SrgbFromRec709)
 637 {
 638   const auto srgb = ColorProfileDesc::From({
 639       Chromaticities::Srgb(),
 640       PiecewiseGammaDesc::Srgb(),
 641   });
 642   const auto rec709 = ColorProfileDesc::From({
 643       Chromaticities::Rec709(),
 644       PiecewiseGammaDesc::Rec709(),
 645   });
 646
 647   {
 648     const auto conv = ColorProfileConversionDesc::From({
 649         .src = srgb,
 650         .dst = srgb,
 651     });
 652     auto src = vec3(16.0);
 653     auto dst = conv.Apply(src / 255) * 255;
 654
 655     const auto tfa = PiecewiseGammaDesc::Srgb();
 656     const auto tfb = PiecewiseGammaDesc::Srgb();
 657     const auto expected =
 658         TfFromLinear(tfb, LinearFromTf(tfa, src.x() / 255)) * 255;
 659
 660     printf("%f %f %f\n", src.x(), src.y(), src.z());
 661     printf("%f %f %f\n", dst.x(), dst.y(), dst.z());
 662     EXPECT_LT(Stats::Diff(dst.data, vec3(expected).data), (Stats::Error{0.42}));
 663   }
 664   {
 665     const auto conv = ColorProfileConversionDesc::From({
 666         .src = rec709,
 667         .dst = rec709,
 668     });
 669     auto src = vec3(16.0);
 670     auto dst = conv.Apply(src / 255) * 255;
 671
 672     const auto tfa = PiecewiseGammaDesc::Rec709();
 673     const auto tfb = PiecewiseGammaDesc::Rec709();
 674     const auto expected =
 675         TfFromLinear(tfb, LinearFromTf(tfa, src.x() / 255)) * 255;
 676
 677     printf("%f %f %f\n", src.x(), src.y(), src.z());
 678     printf("%f %f %f\n", dst.x(), dst.y(), dst.z());
 679     EXPECT_LT(Stats::Diff(dst.data, vec3(expected).data), (Stats::Error{1e-6}));
 680   }
 681   {
 682     const auto conv = ColorProfileConversionDesc::From({
 683         .src = rec709,
 684         .dst = srgb,
 685     });
 686     auto src = vec3(16.0);
 687     auto dst = conv.Apply(src / 255) * 255;
 688
 689     const auto tfa = PiecewiseGammaDesc::Rec709();
 690     const auto tfb = PiecewiseGammaDesc::Srgb();
 691     const auto expected =
 692         TfFromLinear(tfb, LinearFromTf(tfa, src.x() / 255)) * 255;
 693     printf("expected: %f\n", expected);
 694     printf("%f %f %f\n", src.x(), src.y(), src.z());
 695     printf("%f %f %f\n", dst.x(), dst.y(), dst.z());
 696     EXPECT_LT(Stats::Diff(dst.data, vec3(expected).data), (Stats::Error{0.12}));
 697   }
 698 }