gfx/gl/Colorspaces.cpp

   1 /* This Source Code Form is subject to the terms of the Mozilla Public
   2  * License, v. 2.0. If a copy of the MPL was not distributed with this
   3  * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
   4
   5 // We are going to be doing so, so many transforms, so descriptive labels are
   6 // critical.
   7
   8 #include "Colorspaces.h"
   9
  10 #include "nsDebug.h"
  11 #include "qcms.h"
  12
  13 namespace mozilla::color {
  14
  15 // tf = { k * linear                   | linear < b
  16 //      { a * pow(linear, 1/g) - (1-a) | linear >= b
  17 float TfFromLinear(const PiecewiseGammaDesc& desc, const float linear) {
  18   if (linear < desc.b) {
  19     return linear * desc.k;
  20   }
  21   float ret = linear;
  22   ret = powf(ret, 1.0f / desc.g);
  23   ret *= desc.a;
  24   ret -= (desc.a - 1);
  25   return ret;
  26 }
  27
  28 float LinearFromTf(const PiecewiseGammaDesc& desc, const float tf) {
  29   const auto linear_if_low = tf / desc.k;
  30   if (linear_if_low < desc.b) {
  31     return linear_if_low;
  32   }
  33   float ret = tf;
  34   ret += (desc.a - 1);
  35   ret /= desc.a;
  36   ret = powf(ret, 1.0f * desc.g);
  37   return ret;
  38 }
  39
  40 // -
  41
  42 mat3 YuvFromRgb(const YuvLumaCoeffs& yc) {
  43   // Y is always [0,1]
  44   // U and V are signed, and could be either [-1,+1] or [-0.5,+0.5].
  45   // Specs generally use [-0.5,+0.5], so we use that too.
  46   // E.g.
  47   // y = 0.2126*r + 0.7152*g + 0.0722*b
  48   // u = (b - y) / (u_range = u_max - u_min) // u_min = -u_max
  49   //   = (b - y) / (u(0,0,1) - u(1,1,0))
  50   //   = (b - y) / (2 * u(0,0,1))
  51   //   = (b - y) / (2 * u.b))
  52   //   = (b - y) / (2 * (1 - 0.0722))
  53   //   = (-0.2126*r + -0.7152*g + (1-0.0722)*b) / 1.8556
  54   // v = (r - y) / 1.5748;
  55   //   = ((1-0.2126)*r + -0.7152*g + -0.0722*b) / 1.5748
  56   const auto y = vec3({yc.r, yc.g, yc.b});
  57   const auto u = vec3({0, 0, 1}) - y;
  58   const auto v = vec3({1, 0, 0}) - y;
  59
  60   // From rows:
  61   return mat3({y, u / (2 * u.z()), v / (2 * v.x())});
  62 }
  63
  64 mat4 YuvFromYcbcr(const YcbcrDesc& d) {
  65   // E.g.
  66   // y = (yy - 16) / (235 - 16); // 16->0, 235->1
  67   // u = (cb - 128) / (240 - 16); // 16->-0.5, 128->0, 240->+0.5
  68   // v = (cr - 128) / (240 - 16);
  69
  70   const auto yRange = d.y1 - d.y0;
  71   const auto uHalfRange = d.uPlusHalf - d.u0;
  72   const auto uRange = 2 * uHalfRange;
  73
  74   const auto ycbcrFromYuv = mat4{{vec4{{yRange, 0, 0, d.y0}},
  75                                   {{0, uRange, 0, d.u0}},
  76                                   {{0, 0, uRange, d.u0}},
  77                                   {{0, 0, 0, 1}}}};
  78   const auto yuvFromYcbcr = inverse(ycbcrFromYuv);
  79   return yuvFromYcbcr;
  80 }
  81
  82 inline vec3 CIEXYZ_from_CIExyY(const vec2 xy, const float Y = 1) {
  83   const auto xyz = vec3(xy, 1 - xy.x() - xy.y());
  84   const auto XYZ = xyz * (Y / xy.y());
  85   return XYZ;
  86 }
  87
  88 mat3 XyzFromLinearRgb(const Chromaticities& c) {
  89   // http://www.brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html
  90
  91   // Given red (xr, yr), green (xg, yg), blue (xb, yb),
  92   // and whitepoint (XW, YW, ZW)
  93
  94   // [ X ]       [ R ]
  95   // [ Y ] = M x [ G ]
  96   // [ Z ]       [ B ]
  97
  98   //     [ Sr*Xr Sg*Xg Sb*Xb ]
  99   // M = [ Sr*Yr Sg*Yg Sb*Yb ]
 100   //     [ Sr*Zr Sg*Zg Sb*Zb ]
 101
 102   // Xr = xr / yr
 103   // Yr = 1
 104   // Zr = (1 - xr - yr) / yr
 105
 106   // Xg = xg / yg
 107   // Yg = 1
 108   // Zg = (1 - xg - yg) / yg
 109
 110   // Xb = xb / yb
 111   // Yb = 1
 112   // Zb = (1 - xb - yb) / yb
 113
 114   // [ Sr ]   [ Xr Xg Xb ]^-1   [ XW ]
 115   // [ Sg ] = [ Yr Yg Yb ]    x [ YW ]
 116   // [ Sb ]   [ Zr Zg Zb ]      [ ZW ]
 117
 118   const auto xrgb = vec3({c.rx, c.gx, c.bx});
 119   const auto yrgb = vec3({c.ry, c.gy, c.by});
 120
 121   const auto Xrgb = xrgb / yrgb;
 122   const auto Yrgb = vec3(1);
 123   const auto Zrgb = (vec3(1) - xrgb - yrgb) / yrgb;
 124
 125   const auto XYZrgb = mat3({Xrgb, Yrgb, Zrgb});
 126   const auto XYZrgb_inv = inverse(XYZrgb);
 127   const auto XYZwhitepoint = vec3({c.wx, c.wy, 1 - c.wx - c.wy}) / c.wy;
 128   const auto Srgb = XYZrgb_inv * XYZwhitepoint;
 129
 130   const auto M = mat3({Srgb * Xrgb, Srgb * Yrgb, Srgb * Zrgb});
 131   return M;
 132 }
 133
 134 // -
 135 ColorspaceTransform ColorspaceTransform::Create(const ColorspaceDesc& src,
 136                                                 const ColorspaceDesc& dst) {
 137   auto ct = ColorspaceTransform{src, dst};
 138   ct.srcTf = src.tf;
 139   ct.dstTf = dst.tf;
 140
 141   const auto RgbTfFrom = [&](const ColorspaceDesc& cs) {
 142     auto rgbFrom = mat4::Identity();
 143     if (cs.yuv) {
 144       const auto yuvFromYcbcr = YuvFromYcbcr(cs.yuv->ycbcr);
 145       const auto yuvFromRgb = YuvFromRgb(cs.yuv->yCoeffs);
 146       const auto rgbFromYuv = inverse(yuvFromRgb);
 147       const auto rgbFromYuv4 = mat4(rgbFromYuv);
 148
 149       const auto rgbFromYcbcr = rgbFromYuv4 * yuvFromYcbcr;
 150       rgbFrom = rgbFromYcbcr;
 151     }
 152     return rgbFrom;
 153   };
 154
 155   ct.srcRgbTfFromSrc = RgbTfFrom(src);
 156   const auto dstRgbTfFromDst = RgbTfFrom(dst);
 157   ct.dstFromDstRgbTf = inverse(dstRgbTfFromDst);
 158
 159   // -
 160
 161   ct.dstRgbLinFromSrcRgbLin = mat3::Identity();
 162   if (!(src.chrom == dst.chrom)) {
 163     const auto xyzFromSrcRgbLin = XyzFromLinearRgb(src.chrom);
 164     const auto xyzFromDstRgbLin = XyzFromLinearRgb(dst.chrom);
 165     const auto dstRgbLinFromXyz = inverse(xyzFromDstRgbLin);
 166     ct.dstRgbLinFromSrcRgbLin = dstRgbLinFromXyz * xyzFromSrcRgbLin;
 167   }
 168
 169   return ct;
 170 }
 171
 172 vec3 ColorspaceTransform::DstFromSrc(const vec3 src) const {
 173   const auto srcRgbTf = srcRgbTfFromSrc * vec4(src, 1);
 174   auto srcRgbLin = srcRgbTf;
 175   if (srcTf) {
 176     srcRgbLin.x(LinearFromTf(*srcTf, srcRgbTf.x()));
 177     srcRgbLin.y(LinearFromTf(*srcTf, srcRgbTf.y()));
 178     srcRgbLin.z(LinearFromTf(*srcTf, srcRgbTf.z()));
 179   }
 180
 181   const auto dstRgbLin = dstRgbLinFromSrcRgbLin * vec3(srcRgbLin);
 182   auto dstRgbTf = dstRgbLin;
 183   if (dstTf) {
 184     dstRgbTf.x(TfFromLinear(*dstTf, dstRgbLin.x()));
 185     dstRgbTf.y(TfFromLinear(*dstTf, dstRgbLin.y()));
 186     dstRgbTf.z(TfFromLinear(*dstTf, dstRgbLin.z()));
 187   }
 188
 189   const auto dst4 = dstFromDstRgbTf * vec4(dstRgbTf, 1);
 190   return vec3(dst4);
 191 }
 192
 193 // -
 194
 195 mat3 XyzAFromXyzB_BradfordLinear(const vec2 xyA, const vec2 xyB) {
 196   // This is what ICC profiles use to do whitepoint transforms,
 197   // because ICC also requires D50 for the Profile Connection Space.
 198
 199   // From https://www.color.org/specification/ICC.1-2022-05.pdf
 200   // E.3 "Linearized Bradford transformation":
 201
 202   constexpr auto M_BFD = mat3{{
 203       vec3{{0.8951, 0.2664f, -0.1614f}},
 204       vec3{{-0.7502f, 1.7135f, 0.0367f}},
 205       vec3{{0.0389f, -0.0685f, 1.0296f}},
 206   }};
 207   // NB: They use rho/gamma/beta, but we'll use R/G/B here.
 208   const auto XYZDst = CIEXYZ_from_CIExyY(xyA);  // "XYZ_W", WP of PCS
 209   const auto XYZSrc = CIEXYZ_from_CIExyY(xyB);  // "XYZ_NAW", WP of src
 210   const auto rgbSrc = M_BFD * XYZSrc;           // "RGB_SRC"
 211   const auto rgbDst = M_BFD * XYZDst;           // "RGB_PCS"
 212   const auto rgbDstOverSrc = rgbDst / rgbSrc;
 213   const auto M_dstOverSrc = mat3::Scale(rgbDstOverSrc);
 214   const auto M_adapt = inverse(M_BFD) * M_dstOverSrc * M_BFD;
 215   return M_adapt;
 216 }
 217
 218 std::optional<mat4> ColorspaceTransform::ToMat4() const {
 219   mat4 fromSrc = srcRgbTfFromSrc;
 220   if (srcTf) return {};
 221   fromSrc = mat4(dstRgbLinFromSrcRgbLin) * fromSrc;
 222   if (dstTf) return {};
 223   fromSrc = dstFromDstRgbTf * fromSrc;
 224   return fromSrc;
 225 }
 226
 227 Lut3 ColorspaceTransform::ToLut3(const ivec3 size) const {
 228   auto lut = Lut3::Create(size);
 229   lut.SetMap([&](const vec3& srcVal) { return DstFromSrc(srcVal); });
 230   return lut;
 231 }
 232
 233 vec3 Lut3::Sample(const vec3 in01) const {
 234   const auto coord = vec3(size - 1) * in01;
 235   const auto p0 = floor(coord);
 236   const auto dp = coord - p0;
 237   const auto ip0 = ivec3(p0);
 238
 239   // Trilinear
 240   const auto f000 = Fetch(ip0 + ivec3({0, 0, 0}));
 241   const auto f100 = Fetch(ip0 + ivec3({1, 0, 0}));
 242   const auto f010 = Fetch(ip0 + ivec3({0, 1, 0}));
 243   const auto f110 = Fetch(ip0 + ivec3({1, 1, 0}));
 244   const auto f001 = Fetch(ip0 + ivec3({0, 0, 1}));
 245   const auto f101 = Fetch(ip0 + ivec3({1, 0, 1}));
 246   const auto f011 = Fetch(ip0 + ivec3({0, 1, 1}));
 247   const auto f111 = Fetch(ip0 + ivec3({1, 1, 1}));
 248
 249   const auto fx00 = mix(f000, f100, dp.x());
 250   const auto fx10 = mix(f010, f110, dp.x());
 251   const auto fx01 = mix(f001, f101, dp.x());
 252   const auto fx11 = mix(f011, f111, dp.x());
 253
 254   const auto fxy0 = mix(fx00, fx10, dp.y());
 255   const auto fxy1 = mix(fx01, fx11, dp.y());
 256
 257   const auto fxyz = mix(fxy0, fxy1, dp.z());
 258   return fxyz;
 259 }
 260
 261 // -
 262
 263 ColorProfileDesc ColorProfileDesc::From(const ColorspaceDesc& cspace) {
 264   auto ret = ColorProfileDesc{};
 265
 266   if (cspace.yuv) {
 267     const auto yuvFromYcbcr = YuvFromYcbcr(cspace.yuv->ycbcr);
 268     const auto yuvFromRgb = YuvFromRgb(cspace.yuv->yCoeffs);
 269     const auto rgbFromYuv = inverse(yuvFromRgb);
 270     ret.rgbFromYcbcr = mat4(rgbFromYuv) * yuvFromYcbcr;
 271   }
 272
 273   if (cspace.tf) {
 274     const size_t tableSize = 256;
 275     auto& tableR = ret.linearFromTf.r;
 276     tableR.resize(tableSize);
 277     for (size_t i = 0; i < tableR.size(); i++) {
 278       const float tfVal = i / float(tableR.size() - 1);
 279       const float linearVal = LinearFromTf(*cspace.tf, tfVal);
 280       tableR[i] = linearVal;
 281     }
 282     ret.linearFromTf.g = tableR;
 283     ret.linearFromTf.b = tableR;
 284   }
 285
 286   ret.xyzd65FromLinearRgb = XyzFromLinearRgb(cspace.chrom);
 287
 288   return ret;
 289 }
 290
 291 // -
 292
 293 template <class T>
 294 constexpr inline T NewtonEstimateX(const T x1, const T y1, const T dydx,
 295                                    const T y2 = 0) {
 296   // Estimate x s.t. y=0
 297   // y = y0 + x*dydx;
 298   // y0 = y - x*dydx;
 299   // y1 - x1*dydx = y2 - x2*dydx
 300   // x2*dydx = y2 - y1 + x1*dydx
 301   // x2 = (y2 - y1)/dydx + x1
 302   return (y2 - y1) / dydx + x1;
 303 }
 304
 305 float GuessGamma(const std::vector<float>& vals, float exp_guess) {
 306   // Approximate (signed) error = 0.0.
 307   constexpr float d_exp = 0.001;
 308   constexpr float error_tolerance = 0.001;
 309   struct Samples {
 310     float y1, y2;
 311   };
 312   const auto Sample = [&](const float exp) {
 313     int i = -1;
 314     auto samples = Samples{};
 315     for (const auto& expected : vals) {
 316       i += 1;
 317       const auto in = i / float(vals.size() - 1);
 318       samples.y1 += powf(in, exp) - expected;
 319       samples.y2 += powf(in, exp + d_exp) - expected;
 320     }
 321     samples.y1 /= vals.size();  // Normalize by val count.
 322     samples.y2 /= vals.size();
 323     return samples;
 324   };
 325   constexpr int MAX_ITERS = 10;
 326   for (int i = 1;; i++) {
 327     const auto err = Sample(exp_guess);
 328     const auto derr = err.y2 - err.y1;
 329     exp_guess = NewtonEstimateX(exp_guess, err.y1, derr / d_exp);
 330     // Check if we were close before, because then this last round of estimation
 331     // should get us pretty much right on it.
 332     if (std::abs(err.y1) < error_tolerance) {
 333       return exp_guess;
 334     }
 335     if (i >= MAX_ITERS) {
 336       printf_stderr("GuessGamma() -> %f after %i iterations (avg err %f)\n",
 337                     exp_guess, i, err.y1);
 338       MOZ_ASSERT(false, "GuessGamma failed.");
 339       return exp_guess;
 340     }
 341   }
 342 }
 343
 344 // -
 345
 346 ColorProfileDesc ColorProfileDesc::From(const qcms_profile& qcmsProfile) {
 347   ColorProfileDesc ret;
 348
 349   qcms_profile_data data = {};
 350   qcms_profile_get_data(&qcmsProfile, &data);
 351
 352   auto xyzd50FromLinearRgb = mat3{};
 353   // X contributions from [R,G,B]
 354   xyzd50FromLinearRgb.at(0, 0) = data.red_colorant_xyzd50[0];
 355   xyzd50FromLinearRgb.at(1, 0) = data.green_colorant_xyzd50[0];
 356   xyzd50FromLinearRgb.at(2, 0) = data.blue_colorant_xyzd50[0];
 357   // Y contributions from [R,G,B]
 358   xyzd50FromLinearRgb.at(0, 1) = data.red_colorant_xyzd50[1];
 359   xyzd50FromLinearRgb.at(1, 1) = data.green_colorant_xyzd50[1];
 360   xyzd50FromLinearRgb.at(2, 1) = data.blue_colorant_xyzd50[1];
 361   // Z contributions from [R,G,B]
 362   xyzd50FromLinearRgb.at(0, 2) = data.red_colorant_xyzd50[2];
 363   xyzd50FromLinearRgb.at(1, 2) = data.green_colorant_xyzd50[2];
 364   xyzd50FromLinearRgb.at(2, 2) = data.blue_colorant_xyzd50[2];
 365
 366   const auto d65FromD50 = XyzAFromXyzB_BradfordLinear(D65, D50);
 367   ret.xyzd65FromLinearRgb = d65FromD50 * xyzd50FromLinearRgb;
 368
 369   // -
 370
 371   const auto Fn = [&](std::vector<float>* const linearFromTf,
 372                       int32_t claimed_samples,
 373                       const qcms_color_channel channel) {
 374     if (claimed_samples == 0) return;  // No tf.
 375
 376     if (claimed_samples == -1) {
 377       claimed_samples = 4096;  // Ask it to generate a bunch.
 378       claimed_samples = 256;   // Ask it to generate a bunch.
 379     }
 380
 381     linearFromTf->resize(AssertedCast<size_t>(claimed_samples));
 382
 383     const auto begin = linearFromTf->data();
 384     qcms_profile_get_lut(&qcmsProfile, channel, begin,
 385                          begin + linearFromTf->size());
 386   };
 387
 388   Fn(&ret.linearFromTf.r, data.linear_from_trc_red_samples,
 389      qcms_color_channel::Red);
 390   Fn(&ret.linearFromTf.b, data.linear_from_trc_blue_samples,
 391      qcms_color_channel::Blue);
 392   Fn(&ret.linearFromTf.g, data.linear_from_trc_green_samples,
 393      qcms_color_channel::Green);
 394
 395   // -
 396
 397   return ret;
 398 }
 399
 400 // -
 401
 402 ColorProfileConversionDesc ColorProfileConversionDesc::From(
 403     const FromDesc& desc) {
 404   const auto dstLinearRgbFromXyzd65 = inverse(desc.dst.xyzd65FromLinearRgb);
 405   auto ret = ColorProfileConversionDesc{
 406       .srcRgbFromSrcYuv = desc.src.rgbFromYcbcr,
 407       .srcLinearFromSrcTf = desc.src.linearFromTf,
 408       .dstLinearFromSrcLinear =
 409           dstLinearRgbFromXyzd65 * desc.src.xyzd65FromLinearRgb,
 410       .dstTfFromDstLinear = {},
 411   };
 412   bool sameTF = true;
 413   sameTF &= desc.src.linearFromTf.r == desc.dst.linearFromTf.r;
 414   sameTF &= desc.src.linearFromTf.g == desc.dst.linearFromTf.g;
 415   sameTF &= desc.src.linearFromTf.b == desc.dst.linearFromTf.b;
 416   if (sameTF) {
 417     ret.srcLinearFromSrcTf = {};
 418     ret.dstTfFromDstLinear = {};
 419   } else {
 420     const auto Invert = [](const std::vector<float>& linearFromTf,
 421                            std::vector<float>* const tfFromLinear) {
 422       const auto size = linearFromTf.size();
 423       MOZ_ASSERT(size != 1);  // Less than two is uninvertable.
 424       if (size < 2) return;
 425       (*tfFromLinear).resize(size);
 426       InvertLut(linearFromTf, &*tfFromLinear);
 427     };
 428     Invert(desc.dst.linearFromTf.r, &ret.dstTfFromDstLinear.r);
 429     Invert(desc.dst.linearFromTf.g, &ret.dstTfFromDstLinear.g);
 430     Invert(desc.dst.linearFromTf.b, &ret.dstTfFromDstLinear.b);
 431   }
 432   return ret;
 433 }
 434
 435 }  // namespace mozilla::color