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
8 #include "Colorspaces.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
;
22 ret
= powf(ret
, 1.0f
/ desc
.g
);
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
) {
36 ret
= powf(ret
, 1.0f
* desc
.g
);
42 mat3
YuvFromRgb(const YuvLumaCoeffs
& yc
) {
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.
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
;
61 return mat3({y
, u
/ (2 * u
.z()), v
/ (2 * v
.x())});
64 mat4
YuvFromYcbcr(const YcbcrDesc
& d
) {
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
}},
78 const auto yuvFromYcbcr
= inverse(ycbcrFromYuv
);
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());
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)
98 // [ Sr*Xr Sg*Xg Sb*Xb ]
99 // M = [ Sr*Yr Sg*Yg Sb*Yb ]
100 // [ Sr*Zr Sg*Zg Sb*Zb ]
104 // Zr = (1 - xr - yr) / yr
108 // Zg = (1 - xg - yg) / yg
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
});
135 ColorspaceTransform
ColorspaceTransform::Create(const ColorspaceDesc
& src
,
136 const ColorspaceDesc
& dst
) {
137 auto ct
= ColorspaceTransform
{src
, dst
};
141 const auto RgbTfFrom
= [&](const ColorspaceDesc
& cs
) {
142 auto rgbFrom
= mat4::Identity();
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
;
155 ct
.srcRgbTfFromSrc
= RgbTfFrom(src
);
156 const auto dstRgbTfFromDst
= RgbTfFrom(dst
);
157 ct
.dstFromDstRgbTf
= inverse(dstRgbTfFromDst
);
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
;
172 vec3
ColorspaceTransform::DstFromSrc(const vec3 src
) const {
173 const auto srcRgbTf
= srcRgbTfFromSrc
* vec4(src
, 1);
174 auto srcRgbLin
= srcRgbTf
;
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
;
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);
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
;
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
;
227 Lut3
ColorspaceTransform::ToLut3(const ivec3 size
) const {
228 auto lut
= Lut3::Create(size
);
229 lut
.SetMap([&](const vec3
& srcVal
) { return DstFromSrc(srcVal
); });
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
);
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());
263 ColorProfileDesc
ColorProfileDesc::From(const ColorspaceDesc
& cspace
) {
264 auto ret
= ColorProfileDesc
{};
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
;
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
);
294 constexpr inline T
NewtonEstimateX(const T x1
, const T y1
, const T dydx
,
296 // Estimate x s.t. y=0
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;
312 const auto Sample
= [&](const float exp
) {
314 auto samples
= Samples
{};
315 for (const auto& expected
: vals
) {
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();
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
) {
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.");
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
;
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
);
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
= {},
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
;
417 ret
.srcLinearFromSrcTf
= {};
418 ret
.dstTfFromDstLinear
= {};
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
);
435 } // namespace mozilla::color