no bug - Correct some typos in the comments. a=typo-fix
[gecko.git] / gfx / gl / Colorspaces.cpp
blob29ddf37cd9a6ac249a3a57a2c134f9566de2ff24
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/. */
5 // We are going to be doing so, so many transforms, so descriptive labels are
6 // critical.
8 #include "Colorspaces.h"
10 #include "nsDebug.h"
11 #include "qcms.h"
13 namespace mozilla::color {
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;
21 float ret = linear;
22 ret = powf(ret, 1.0f / desc.g);
23 ret *= desc.a;
24 ret -= (desc.a - 1);
25 return ret;
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;
33 float ret = tf;
34 ret += (desc.a - 1);
35 ret /= desc.a;
36 ret = powf(ret, 1.0f * desc.g);
37 return ret;
40 // -
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;
60 // From rows:
61 return mat3({y, u / (2 * u.z()), v / (2 * v.x())});
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);
70 const auto yRange = d.y1 - d.y0;
71 const auto uHalfRange = d.uPlusHalf - d.u0;
72 const auto uRange = 2 * uHalfRange;
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;
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;
88 mat3 XyzFromLinearRgb(const Chromaticities& c) {
89 // http://www.brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html
91 // Given red (xr, yr), green (xg, yg), blue (xb, yb),
92 // and whitepoint (XW, YW, ZW)
94 // [ X ] [ R ]
95 // [ Y ] = M x [ G ]
96 // [ Z ] [ B ]
98 // [ Sr*Xr Sg*Xg Sb*Xb ]
99 // M = [ Sr*Yr Sg*Yg Sb*Yb ]
100 // [ Sr*Zr Sg*Zg Sb*Zb ]
102 // Xr = xr / yr
103 // Yr = 1
104 // Zr = (1 - xr - yr) / yr
106 // Xg = xg / yg
107 // Yg = 1
108 // Zg = (1 - xg - yg) / yg
110 // Xb = xb / yb
111 // Yb = 1
112 // Zb = (1 - xb - yb) / yb
114 // [ Sr ] [ Xr Xg Xb ]^-1 [ XW ]
115 // [ Sg ] = [ Yr Yg Yb ] x [ YW ]
116 // [ Sb ] [ Zr Zg Zb ] [ ZW ]
118 const auto xrgb = vec3({c.rx, c.gx, c.bx});
119 const auto yrgb = vec3({c.ry, c.gy, c.by});
121 const auto Xrgb = xrgb / yrgb;
122 const auto Yrgb = vec3(1);
123 const auto Zrgb = (vec3(1) - xrgb - yrgb) / yrgb;
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;
130 const auto M = mat3({Srgb * Xrgb, Srgb * Yrgb, Srgb * Zrgb});
131 return M;
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;
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);
149 const auto rgbFromYcbcr = rgbFromYuv4 * yuvFromYcbcr;
150 rgbFrom = rgbFromYcbcr;
152 return rgbFrom;
155 ct.srcRgbTfFromSrc = RgbTfFrom(src);
156 const auto dstRgbTfFromDst = RgbTfFrom(dst);
157 ct.dstFromDstRgbTf = inverse(dstRgbTfFromDst);
159 // -
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;
169 return ct;
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()));
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()));
189 const auto dst4 = dstFromDstRgbTf * vec4(dstRgbTf, 1);
190 return vec3(dst4);
193 // -
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.
199 // From https://www.color.org/specification/ICC.1-2022-05.pdf
200 // E.3 "Linearized Bradford transformation":
202 const 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}},
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;
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;
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;
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);
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}));
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());
254 const auto fxy0 = mix(fx00, fx10, dp.y());
255 const auto fxy1 = mix(fx01, fx11, dp.y());
257 const auto fxyz = mix(fxy0, fxy1, dp.z());
258 return fxyz;
261 // -
263 ColorProfileDesc ColorProfileDesc::From(const ColorspaceDesc& cspace) {
264 auto ret = ColorProfileDesc{};
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;
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;
282 ret.linearFromTf.g = tableR;
283 ret.linearFromTf.b = tableR;
286 ret.xyzd65FromLinearRgb = XyzFromLinearRgb(cspace.chrom);
288 return ret;
291 // -
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;
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;
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;
321 samples.y1 /= vals.size(); // Normalize by val count.
322 samples.y2 /= vals.size();
323 return samples;
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;
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;
344 // -
346 ColorProfileDesc ColorProfileDesc::From(const qcms_profile& qcmsProfile) {
347 ColorProfileDesc ret;
349 qcms_profile_data data = {};
350 qcms_profile_get_data(&qcmsProfile, &data);
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];
366 const auto d65FromD50 = XyzAFromXyzB_BradfordLinear(D65, D50);
367 ret.xyzd65FromLinearRgb = d65FromD50 * xyzd50FromLinearRgb;
369 // -
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.
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.
381 linearFromTf->resize(AssertedCast<size_t>(claimed_samples));
383 const auto begin = linearFromTf->data();
384 qcms_profile_get_lut(&qcmsProfile, channel, begin,
385 begin + linearFromTf->size());
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);
395 // -
397 return ret;
400 // -
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 = {},
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);
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);
432 return ret;
435 } // namespace mozilla::color