1 /*
2 * Copyright (c) 2015 Kevin Wheatley <kevin.j.wheatley@gmail.com>
3 * Copyright (c) 2016 Ronald S. Bultje <rsbultje@gmail.com>
4 * Copyright (c) 2023 Leo Izen <leo.izen@gmail.com>
5 *
6 * This file is part of FFmpeg.
7 *
8 * FFmpeg is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU Lesser General Public
10 * License as published by the Free Software Foundation; either
11 * version 2.1 of the License, or (at your option) any later version.
12 *
13 * FFmpeg is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 * Lesser General Public License for more details.
17 *
18 * You should have received a copy of the GNU Lesser General Public
19 * License along with FFmpeg; if not, write to the Free Software
20 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
21 */
22
23 /**
24 * @file Colorspace functions for libavutil
25 * @author Ronald S. Bultje <rsbultje@gmail.com>
26 * @author Leo Izen <leo.izen@gmail.com>
27 * @author Kevin Wheatley <kevin.j.wheatley@gmail.com>
28 */
29
30 #include <stdlib.h>
31 #include <math.h>
32
37
38 #define AVR(d) { (int)(d * 100000 + 0.5), 100000 }
39
40 /*
41 * All constants explained in e.g. https://linuxtv.org/downloads/v4l-dvb-apis/ch02s06.html
42 * The older ones (bt470bg/m) are also explained in their respective ITU docs
43 * (e.g. https://www.itu.int/dms_pubrec/itu-r/rec/bt/R-REC-BT.470-5-199802-S!!PDF-E.pdf)
44 * whereas the newer ones can typically be copied directly from wikipedia :)
45 */
56 };
57
59 {
61
67
68 return coeffs;
69 }
70
71 #define WP_D65 { AVR(0.3127), AVR(0.3290) }
72 #define WP_C { AVR(0.3100), AVR(0.3160) }
73 #define WP_DCI { AVR(0.3140), AVR(0.3510) }
74 #define WP_E { {1, 3}, {1, 3} }
75
88 };
89
91 {
93
99
101 }
102
104 {
106 /* denominator assumed to be positive */
108 }
109
111 {
113
116 if (!
ref->prim.r.x.num)
117 continue;
118
127
130 }
131
133 }
134
147 };
148
150 {
151 double gamma;
153 return 0.0;
155 if (gamma > 0)
156 return gamma;
157 return 0.0;
158 }
159
160 #define BT709_alpha 1.099296826809442
161 #define BT709_beta 0.018053968510807
162
164 {
167
168 return (0.0 > Lc) ? 0.0
169 : (
b > Lc) ? 4.500 * Lc
170 :
a * pow(Lc, 0.45) - (
a - 1.0);
171 }
172
174 {
177
178 return (0.0 >
E) ? 0.0
179 : (
b >
E) ?
E / 4.500
180 : pow((
E + (
a - 1.0)) /
a, 1.0 / 0.45);
181 }
182
184 {
185 return (0.0 > Lc) ? 0.0 : pow(Lc, 1.0/ 2.2);
186 }
187
189 {
190 return (0.0 >
E) ? 0.0 : pow(
E, 2.2);
191 }
192
194 {
195 return (0.0 > Lc) ? 0.0 : pow(Lc, 1.0/ 2.8);
196 }
197
199 {
200 return (0.0 >
E) ? 0.0 : pow(
E, 2.8);
201 }
202
204 {
205 const double a = 1.1115;
206 const double b = 0.0228;
207
208 return (0.0 > Lc) ? 0.0
209 : (
b > Lc) ? 4.000 * Lc
210 :
a * pow(Lc, 0.45) - (
a - 1.0);
211 }
212
214 {
215 const double a = 1.1115;
216 const double b = 4.000 * 0.0228;
217
218 return (0.0 >
E) ? 0.0
219 : (
b >
E) ?
E / 4.000
220 : pow((
E + (
a - 1.0)) /
a, 1.0 / 0.45);
221 }
222
224 {
225 return Lc;
226 }
227
229 {
230 return (0.01 > Lc) ? 0.0 : 1.0 + log10(Lc) / 2.0;
231 }
232
234 {
235 return (0.0 >
E) ? 0.01 : pow(10.0, 2.0 * (
E - 1.0));
236 }
237
239 {
240 // sqrt(10) / 1000
241 return (0.00316227766 > Lc) ? 0.0 : 1.0 + log10(Lc) / 2.5;
242 }
243
245 {
246 return (0.0 >
E) ? 0.00316227766 : pow(10.0, 2.5 * (
E - 1.0));
247 }
248
250 {
253
254 return (-
b >= Lc) ? -
a * pow(-Lc, 0.45) + (
a - 1.0)
255 : (
b > Lc) ? 4.500 * Lc
256 :
a * pow( Lc, 0.45) - (
a - 1.0);
257 }
258
260 {
263
264 return (-
b >=
E) ? -pow((-
E + (
a - 1.0)) /
a, 1.0 / 0.45)
265 : (
b >
E) ?
E / 4.500
266 : pow((
E + (
a - 1.0)) /
a, 1.0 / 0.45);
267 }
268
270 {
273
274 return (-0.0045 >= Lc) ? -(
a * pow(-4.0 * Lc, 0.45) + (
a - 1.0)) / 4.0
275 : (
b > Lc) ? 4.500 * Lc
276 :
a * pow( Lc, 0.45) - (
a - 1.0);
277 }
278
280 {
283
284 return (-0.02025 >=
E) ? -pow((-4.0 *
E - (
a - 1.0)) /
a, 1.0 / 0.45) / 4.0
285 : (
b >
E) ?
E / 4.500
286 : pow((
E + (
a - 1.0)) /
a, 1.0 / 0.45);
287 }
288
290 {
291 const double a = 1.055;
292 const double b = 0.0031308;
293
294 return (0.0 > Lc) ? 0.0
295 : (
b > Lc) ? 12.92 * Lc
296 :
a * pow(Lc, 1.0 / 2.4) - (
a - 1.0);
297 }
298
300 {
301 const double a = 1.055;
302 const double b = 12.92 * 0.0031308;
303
304 return (0.0 >
E) ? 0.0
305 : (
b >
E) ?
E / 12.92
306 : pow((
E + (
a - 1.0)) /
a, 2.4);
308 }
309
310 #define PQ_c1 ( 3424.0 / 4096.0) /* c3-c2 + 1 */
311 #define PQ_c2 ( 32.0 * 2413.0 / 4096.0)
312 #define PQ_c3 ( 32.0 * 2392.0 / 4096.0)
313 #define PQ_m (128.0 * 2523.0 / 4096.0)
314 #define PQ_n ( 0.25 * 2610.0 / 4096.0)
315
317 {
320 const double c3 =
PQ_c3;
321 const double m =
PQ_m;
322 const double n =
PQ_n;
323 const double L = Lc / 10000.0;
324 const double Ln = pow(
L, n);
325
326 return (0.0 > Lc) ? 0.0
327 : pow((
c1 +
c2 * Ln) / (1.0 + c3 * Ln), m);
328
329 }
330
332 {
335 const double c3 =
PQ_c3;
336 const double m =
PQ_m;
337 const double n =
PQ_n;
338 const double Em = pow(
E, 1.0 / m);
339
340 return (
c1 > Em) ? 0.0
341 : 10000.0 * pow((Em -
c1) / (
c2 - c3 * Em), 1.0 / n);
342 }
343
344 #define DCI_L 48.00
345 #define DCI_P 52.37
346
348 {
349 return (0.0 > Lc) ? 0.0 : pow(
DCI_L /
DCI_P * Lc, 1.0 / 2.6);
350 }
351
353 {
355 }
356
357 #define HLG_a 0.17883277
358 #define HLG_b 0.28466892
359 #define HLG_c 0.55991073
360
362 // The function uses the definition from HEVC, which assumes that the peak
363 // white is input level = 1. (this is equivalent to scaling E = Lc * 12 and
364 // using the definition from the ARIB STD-B67 spec)
368 return (0.0 > Lc) ? 0.0 :
369 (Lc <= 1.0 / 12.0 ? sqrt(3.0 * Lc) :
a *
log(12.0 * Lc -
b) +
c);
370 }
371
373 {
377 return (0.0 >
E) ? 0.0 :
378 (
E <= 0.5 ?
E *
E / 3.0 : (
exp((
E -
c) /
a) +
b) / 12.0);
379 }
380
398 };
399
401 {
405 }
406
424 };
425
427 {
431 }
432
434 {
435 for (
int i = 0;
i < 3;
i++)
436 E[
i] = (Lw - Lb) *
E[
i] + Lb;
437 }
438
440 {
441 for (
int i = 0;
i < 3;
i++)
442 L[
i] = (
L[
i] - Lb) / (Lw - Lb);
443 }
444
445 #define WRAP_SDR_OETF(name) \
446 static void oetf_##name(double L[3]) \
447 { \
448 for (int i = 0; i < 3; i++) \
449 L[i] = trc_##name(L[i]); \
450 } \
451 \
452 static void oetf_##name##_inv(double E[3]) \
453 { \
454 for (int i = 0; i < 3; i++) \
455 E[i] = trc_##name##_inv(E[i]); \
456 }
457
461
462 #define WRAP_SDR_EOTF(name) \
463 static void eotf_##name(double Lw, double Lb, double E[3]) \
464 { \
465 oetf_##name##_inv(E); \
466 eotf_linear(Lw, Lb, E); \
467 } \
468 \
469 static void eotf_##name##_inv(double Lw, double Lb, double L[3]) \
470 { \
471 eotf_linear_inv(Lw, Lb, L); \
472 oetf_##name(L); \
473 }
474
478
480 {
481 const double Lw_inv = pow(Lw, 1.0 / 2.4);
482 const double Lb_inv = pow(Lb, 1.0 / 2.4);
483 const double a = pow(Lw_inv - Lb_inv, 2.4);
484 const double b = Lb_inv / (Lw_inv - Lb_inv);
485
486 for (
int i = 0;
i < 3;
i++)
487 E[
i] = (-
b >
E[
i]) ? 0.0 :
a * pow(
E[
i] +
b, 2.4);
488 }
489
491 {
492 const double Lw_inv = pow(Lw, 1.0 / 2.4);
493 const double Lb_inv = pow(Lb, 1.0 / 2.4);
494 const double a = pow(Lw_inv - Lb_inv, 2.4);
495 const double b = Lb_inv / (Lw_inv - Lb_inv);
496
497 for (
int i = 0;
i < 3;
i++)
498 L[
i] = (0.0 >
L[
i]) ? 0.0 : pow(
L[
i] /
a, 1.0 / 2.4) -
b;
499 }
500
502 {
503 for (
int i = 0;
i < 3;
i++)
505 }
506
508 {
509 for (
int i = 0;
i < 3;
i++)
511 }
512
513 /* This implementation assumes an SMPTE RP 431-2 reference projector (DCI) */
516 #define DCI_X (42.94 / DCI_L)
517 #define DCI_Z (45.82 / DCI_L)
518
520 {
521 const double Lw[3] = {
DCI_X * Lw_Y, Lw_Y,
DCI_Z * Lw_Y };
522 const double Lb[3] = {
DCI_X * Lb_Y, Lb_Y,
DCI_Z * Lb_Y };
523
524 for (
int i = 0;
i < 3;
i++) {
526 E[
i] =
E[
i] * (Lw[
i] - Lb[
i]) + Lb[
i];
527 }
528 }
529
531 {
532 const double Lw[3] = {
DCI_X * Lw_Y, Lw_Y,
DCI_Z * Lw_Y };
533 const double Lb[3] = {
DCI_X * Lb_Y, Lb_Y,
DCI_Z * Lb_Y };
534
535 for (
int i = 0;
i < 3;
i++) {
536 L[
i] = (
L[
i] - Lb[
i]) / (Lw[
i] - Lb[
i]);
538 }
539 }
540
542 {
543 const double gamma =
fmax(1.2 + 0.42 * log10(Lw / 1000.0), 1.0);
544
545 /**
546 * Note: This equation is technically only accurate if the contrast ratio
547 * Lw:Lb is greater than 12:1; otherwise we would need to use a different,
548 * significantly more complicated solution. Ignore this as a highly
549 * degenerate case, since any real world reference display will have a
550 * static contrast ratio multiple orders of magnitude higher.
551 */
552 const double beta = sqrt(3 * pow(Lb / Lw, 1.0 / gamma));
553 double luma;
554
555 for (
int i = 0;
i < 3;
i++)
557
558 luma = 0.2627 *
E[0] + 0.6780 *
E[1] + 0.0593 *
E[2];
559 luma = pow(
fmax(luma, 0.0), gamma - 1.0);
560 for (
int i = 0;
i < 3;
i++)
562 }
563
565 {
566 const double gamma =
fmax(1.2 + 0.42 * log10(Lw / 1000.0), 1.0);
567 const double beta = sqrt(3 * pow(Lb / Lw, 1 / gamma));
568 double luma = 0.2627 *
L[0] + 0.6780 *
L[1] + 0.0593 *
L[2];
569
570 if (luma > 0.0) {
571 luma = pow(luma / Lw, (1 - gamma) / gamma);
572 for (
int i = 0;
i < 3;
i++)
574 } else {
575 L[0] =
L[1] =
L[2] = 0.0;
576 }
577
578 for (
int i = 0;
i < 3;
i++)
580 }
581
589 /* There is no EOTF associated with these logarithmic encodings, since they
590 * are defined purely for transmission of scene referred data. */
593 /* BT.1886 is already defined for values below 0.0, as far as physically
594 * meaningful, so we can directly use it for extended range encodings */
603 };
604
606 {
610 }
611
629 };
630
632 {
636 }