FFmpeg
palette.c
Go to the documentation of this file.
1 /*
2  * Copyright (c) 2020 Björn Ottosson
3  * Copyright (c) 2022 Clément Bœsch <u pkh me>
4  *
5  * This file is part of FFmpeg.
6  *
7  * FFmpeg is free software; you can redistribute it and/or
8  * modify it under the terms of the GNU Lesser General Public
9  * License as published by the Free Software Foundation; either
10  * version 2.1 of the License, or (at your option) any later version.
11  *
12  * FFmpeg is distributed in the hope that it will be useful,
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15  * Lesser General Public License for more details.
16  *
17  * You should have received a copy of the GNU Lesser General Public
18  * License along with FFmpeg; if not, write to the Free Software
19  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
20  */
21 
22 #include "libavutil/common.h"
23 #include "palette.h"
24 
25 #define K ((1 << 16) - 1)
26 #define K2 ((int64_t)K*K)
27 #define P ((1 << 9) - 1)
28 
29 /**
30  * Table mapping formula:
31  * f(x) = x < 0.04045 ? x/12.92 : ((x+0.055)/1.055)^2.4 (sRGB EOTF)
32  * Where x is the normalized index in the table and f(x) the value in the table.
33  * f(x) is remapped to [0;K] and rounded.
34  */
35 static const uint16_t srgb2linear[256] = {
36  0x0000, 0x0014, 0x0028, 0x003c, 0x0050, 0x0063, 0x0077, 0x008b,
37  0x009f, 0x00b3, 0x00c7, 0x00db, 0x00f1, 0x0108, 0x0120, 0x0139,
38  0x0154, 0x016f, 0x018c, 0x01ab, 0x01ca, 0x01eb, 0x020e, 0x0232,
39  0x0257, 0x027d, 0x02a5, 0x02ce, 0x02f9, 0x0325, 0x0353, 0x0382,
40  0x03b3, 0x03e5, 0x0418, 0x044d, 0x0484, 0x04bc, 0x04f6, 0x0532,
41  0x056f, 0x05ad, 0x05ed, 0x062f, 0x0673, 0x06b8, 0x06fe, 0x0747,
42  0x0791, 0x07dd, 0x082a, 0x087a, 0x08ca, 0x091d, 0x0972, 0x09c8,
43  0x0a20, 0x0a79, 0x0ad5, 0x0b32, 0x0b91, 0x0bf2, 0x0c55, 0x0cba,
44  0x0d20, 0x0d88, 0x0df2, 0x0e5e, 0x0ecc, 0x0f3c, 0x0fae, 0x1021,
45  0x1097, 0x110e, 0x1188, 0x1203, 0x1280, 0x1300, 0x1381, 0x1404,
46  0x1489, 0x1510, 0x159a, 0x1625, 0x16b2, 0x1741, 0x17d3, 0x1866,
47  0x18fb, 0x1993, 0x1a2c, 0x1ac8, 0x1b66, 0x1c06, 0x1ca7, 0x1d4c,
48  0x1df2, 0x1e9a, 0x1f44, 0x1ff1, 0x20a0, 0x2150, 0x2204, 0x22b9,
49  0x2370, 0x242a, 0x24e5, 0x25a3, 0x2664, 0x2726, 0x27eb, 0x28b1,
50  0x297b, 0x2a46, 0x2b14, 0x2be3, 0x2cb6, 0x2d8a, 0x2e61, 0x2f3a,
51  0x3015, 0x30f2, 0x31d2, 0x32b4, 0x3399, 0x3480, 0x3569, 0x3655,
52  0x3742, 0x3833, 0x3925, 0x3a1a, 0x3b12, 0x3c0b, 0x3d07, 0x3e06,
53  0x3f07, 0x400a, 0x4110, 0x4218, 0x4323, 0x4430, 0x453f, 0x4651,
54  0x4765, 0x487c, 0x4995, 0x4ab1, 0x4bcf, 0x4cf0, 0x4e13, 0x4f39,
55  0x5061, 0x518c, 0x52b9, 0x53e9, 0x551b, 0x5650, 0x5787, 0x58c1,
56  0x59fe, 0x5b3d, 0x5c7e, 0x5dc2, 0x5f09, 0x6052, 0x619e, 0x62ed,
57  0x643e, 0x6591, 0x66e8, 0x6840, 0x699c, 0x6afa, 0x6c5b, 0x6dbe,
58  0x6f24, 0x708d, 0x71f8, 0x7366, 0x74d7, 0x764a, 0x77c0, 0x7939,
59  0x7ab4, 0x7c32, 0x7db3, 0x7f37, 0x80bd, 0x8246, 0x83d1, 0x855f,
60  0x86f0, 0x8884, 0x8a1b, 0x8bb4, 0x8d50, 0x8eef, 0x9090, 0x9235,
61  0x93dc, 0x9586, 0x9732, 0x98e2, 0x9a94, 0x9c49, 0x9e01, 0x9fbb,
62  0xa179, 0xa339, 0xa4fc, 0xa6c2, 0xa88b, 0xaa56, 0xac25, 0xadf6,
63  0xafca, 0xb1a1, 0xb37b, 0xb557, 0xb737, 0xb919, 0xbaff, 0xbce7,
64  0xbed2, 0xc0c0, 0xc2b1, 0xc4a5, 0xc69c, 0xc895, 0xca92, 0xcc91,
65  0xce94, 0xd099, 0xd2a1, 0xd4ad, 0xd6bb, 0xd8cc, 0xdae0, 0xdcf7,
66  0xdf11, 0xe12e, 0xe34e, 0xe571, 0xe797, 0xe9c0, 0xebec, 0xee1b,
67  0xf04d, 0xf282, 0xf4ba, 0xf6f5, 0xf933, 0xfb74, 0xfdb8, 0xffff,
68 };
69 
70 /**
71  * Table mapping formula:
72  * f(x) = x < 0.0031308 ? x*12.92 : 1.055*x^(1/2.4)-0.055 (sRGB OETF)
73  * Where x is the normalized index in the table and f(x) the value in the table.
74  * f(x) is remapped to [0;0xff] and rounded.
75  *
76  * Since a 16-bit table is too large, we reduce its precision to 9-bit.
77  */
78 static const uint8_t linear2srgb[P + 1] = {
79  0x00, 0x06, 0x0d, 0x12, 0x16, 0x19, 0x1c, 0x1f, 0x22, 0x24, 0x26, 0x28, 0x2a, 0x2c, 0x2e, 0x30,
80  0x32, 0x33, 0x35, 0x36, 0x38, 0x39, 0x3b, 0x3c, 0x3d, 0x3e, 0x40, 0x41, 0x42, 0x43, 0x45, 0x46,
81  0x47, 0x48, 0x49, 0x4a, 0x4b, 0x4c, 0x4d, 0x4e, 0x4f, 0x50, 0x51, 0x52, 0x53, 0x54, 0x55, 0x56,
82  0x56, 0x57, 0x58, 0x59, 0x5a, 0x5b, 0x5b, 0x5c, 0x5d, 0x5e, 0x5f, 0x5f, 0x60, 0x61, 0x62, 0x62,
83  0x63, 0x64, 0x65, 0x65, 0x66, 0x67, 0x67, 0x68, 0x69, 0x6a, 0x6a, 0x6b, 0x6c, 0x6c, 0x6d, 0x6e,
84  0x6e, 0x6f, 0x6f, 0x70, 0x71, 0x71, 0x72, 0x73, 0x73, 0x74, 0x74, 0x75, 0x76, 0x76, 0x77, 0x77,
85  0x78, 0x79, 0x79, 0x7a, 0x7a, 0x7b, 0x7b, 0x7c, 0x7d, 0x7d, 0x7e, 0x7e, 0x7f, 0x7f, 0x80, 0x80,
86  0x81, 0x81, 0x82, 0x82, 0x83, 0x84, 0x84, 0x85, 0x85, 0x86, 0x86, 0x87, 0x87, 0x88, 0x88, 0x89,
87  0x89, 0x8a, 0x8a, 0x8b, 0x8b, 0x8c, 0x8c, 0x8c, 0x8d, 0x8d, 0x8e, 0x8e, 0x8f, 0x8f, 0x90, 0x90,
88  0x91, 0x91, 0x92, 0x92, 0x93, 0x93, 0x93, 0x94, 0x94, 0x95, 0x95, 0x96, 0x96, 0x97, 0x97, 0x97,
89  0x98, 0x98, 0x99, 0x99, 0x9a, 0x9a, 0x9a, 0x9b, 0x9b, 0x9c, 0x9c, 0x9c, 0x9d, 0x9d, 0x9e, 0x9e,
90  0x9f, 0x9f, 0x9f, 0xa0, 0xa0, 0xa1, 0xa1, 0xa1, 0xa2, 0xa2, 0xa3, 0xa3, 0xa3, 0xa4, 0xa4, 0xa5,
91  0xa5, 0xa5, 0xa6, 0xa6, 0xa6, 0xa7, 0xa7, 0xa8, 0xa8, 0xa8, 0xa9, 0xa9, 0xa9, 0xaa, 0xaa, 0xab,
92  0xab, 0xab, 0xac, 0xac, 0xac, 0xad, 0xad, 0xae, 0xae, 0xae, 0xaf, 0xaf, 0xaf, 0xb0, 0xb0, 0xb0,
93  0xb1, 0xb1, 0xb1, 0xb2, 0xb2, 0xb3, 0xb3, 0xb3, 0xb4, 0xb4, 0xb4, 0xb5, 0xb5, 0xb5, 0xb6, 0xb6,
94  0xb6, 0xb7, 0xb7, 0xb7, 0xb8, 0xb8, 0xb8, 0xb9, 0xb9, 0xb9, 0xba, 0xba, 0xba, 0xbb, 0xbb, 0xbb,
95  0xbc, 0xbc, 0xbc, 0xbd, 0xbd, 0xbd, 0xbe, 0xbe, 0xbe, 0xbf, 0xbf, 0xbf, 0xc0, 0xc0, 0xc0, 0xc1,
96  0xc1, 0xc1, 0xc1, 0xc2, 0xc2, 0xc2, 0xc3, 0xc3, 0xc3, 0xc4, 0xc4, 0xc4, 0xc5, 0xc5, 0xc5, 0xc6,
97  0xc6, 0xc6, 0xc6, 0xc7, 0xc7, 0xc7, 0xc8, 0xc8, 0xc8, 0xc9, 0xc9, 0xc9, 0xc9, 0xca, 0xca, 0xca,
98  0xcb, 0xcb, 0xcb, 0xcc, 0xcc, 0xcc, 0xcc, 0xcd, 0xcd, 0xcd, 0xce, 0xce, 0xce, 0xce, 0xcf, 0xcf,
99  0xcf, 0xd0, 0xd0, 0xd0, 0xd0, 0xd1, 0xd1, 0xd1, 0xd2, 0xd2, 0xd2, 0xd2, 0xd3, 0xd3, 0xd3, 0xd4,
100  0xd4, 0xd4, 0xd4, 0xd5, 0xd5, 0xd5, 0xd6, 0xd6, 0xd6, 0xd6, 0xd7, 0xd7, 0xd7, 0xd7, 0xd8, 0xd8,
101  0xd8, 0xd9, 0xd9, 0xd9, 0xd9, 0xda, 0xda, 0xda, 0xda, 0xdb, 0xdb, 0xdb, 0xdc, 0xdc, 0xdc, 0xdc,
102  0xdd, 0xdd, 0xdd, 0xdd, 0xde, 0xde, 0xde, 0xde, 0xdf, 0xdf, 0xdf, 0xe0, 0xe0, 0xe0, 0xe0, 0xe1,
103  0xe1, 0xe1, 0xe1, 0xe2, 0xe2, 0xe2, 0xe2, 0xe3, 0xe3, 0xe3, 0xe3, 0xe4, 0xe4, 0xe4, 0xe4, 0xe5,
104  0xe5, 0xe5, 0xe5, 0xe6, 0xe6, 0xe6, 0xe6, 0xe7, 0xe7, 0xe7, 0xe7, 0xe8, 0xe8, 0xe8, 0xe8, 0xe9,
105  0xe9, 0xe9, 0xe9, 0xea, 0xea, 0xea, 0xea, 0xeb, 0xeb, 0xeb, 0xeb, 0xec, 0xec, 0xec, 0xec, 0xed,
106  0xed, 0xed, 0xed, 0xee, 0xee, 0xee, 0xee, 0xef, 0xef, 0xef, 0xef, 0xef, 0xf0, 0xf0, 0xf0, 0xf0,
107  0xf1, 0xf1, 0xf1, 0xf1, 0xf2, 0xf2, 0xf2, 0xf2, 0xf3, 0xf3, 0xf3, 0xf3, 0xf3, 0xf4, 0xf4, 0xf4,
108  0xf4, 0xf5, 0xf5, 0xf5, 0xf5, 0xf6, 0xf6, 0xf6, 0xf6, 0xf6, 0xf7, 0xf7, 0xf7, 0xf7, 0xf8, 0xf8,
109  0xf8, 0xf8, 0xf9, 0xf9, 0xf9, 0xf9, 0xf9, 0xfa, 0xfa, 0xfa, 0xfa, 0xfb, 0xfb, 0xfb, 0xfb, 0xfb,
110  0xfc, 0xfc, 0xfc, 0xfc, 0xfd, 0xfd, 0xfd, 0xfd, 0xfd, 0xfe, 0xfe, 0xfe, 0xfe, 0xff, 0xff, 0xff,
111 };
112 
113 int32_t ff_srgb_u8_to_linear_int(uint8_t x)
114 {
115  return (int32_t)srgb2linear[x];
116 }
117 
118 uint8_t ff_linear_int_to_srgb_u8(int32_t x)
119 {
120  if (x <= 0) {
121  return 0;
122  } else if (x >= K) {
123  return 0xff;
124  } else {
125  const int32_t xP = x * P;
126  const int32_t i = xP / K;
127  const int32_t m = xP % K;
128  const int32_t y0 = linear2srgb[i];
129  const int32_t y1 = linear2srgb[i + 1];
130  return (m * (y1 - y0) + K/2) / K + y0;
131  }
132 }
133 
134 /* Integer cube root, working only within [0;1] */
135 static int32_t cbrt01_int(int32_t x)
136 {
137  int64_t u;
138 
139  /* Approximation curve is for the [0;1] range */
140  if (x <= 0) return 0;
141  if (x >= K) return K;
142 
143  /*
144  * Initial approximation: x³ - 2.19893x² + 2.01593x + 0.219407
145  *
146  * We are not using any rounding here since the precision is not important
147  * at this stage and it would require the more expensive rounding function
148  * that deals with negative numbers.
149  */
150  u = x*(x*(x + -144107LL) / K + 132114LL) / K + 14379LL;
151 
152  /*
153  * Refine with 2 Halley iterations:
154  * uₙ₊₁ = uₙ-2f(uₙ)f'(uₙ)/(2f'(uₙ)²-f(uₙ)f"(uₙ))
155  * = uₙ(2x+uₙ³)/(x+2uₙ³)
156  *
157  * Note: u is not expected to be < 0, so we can use the (a+b/2)/b rounding.
158  */
159  for (int i = 0; i < 2; i++) {
160  const int64_t u3 = u*u*u;
161  const int64_t den = x + (2*u3 + K2/2) / K2;
162  u = (u * (2*x + (u3 + K2/2) / K2) + den/2) / den;
163  }
164 
165  return u;
166 }
167 
168 static int64_t div_round64(int64_t a, int64_t b) { return (a^b)<0 ? (a-b/2)/b : (a+b/2)/b; }
169 
170 struct Lab ff_srgb_u8_to_oklab_int(uint32_t srgb)
171 {
172  const int32_t r = (int32_t)srgb2linear[srgb >> 16 & 0xff];
173  const int32_t g = (int32_t)srgb2linear[srgb >> 8 & 0xff];
174  const int32_t b = (int32_t)srgb2linear[srgb & 0xff];
175 
176  // Note: lms can actually be slightly over K due to rounded coefficients
177  const int32_t l = (27015LL*r + 35149LL*g + 3372LL*b + K/2) / K;
178  const int32_t m = (13887LL*r + 44610LL*g + 7038LL*b + K/2) / K;
179  const int32_t s = ( 5787LL*r + 18462LL*g + 41286LL*b + K/2) / K;
180 
181  const int32_t l_ = cbrt01_int(l);
182  const int32_t m_ = cbrt01_int(m);
183  const int32_t s_ = cbrt01_int(s);
184 
185  const struct Lab ret = {
186  .L = div_round64( 13792LL*l_ + 52010LL*m_ - 267LL*s_, K),
187  .a = div_round64(129628LL*l_ - 159158LL*m_ + 29530LL*s_, K),
188  .b = div_round64( 1698LL*l_ + 51299LL*m_ - 52997LL*s_, K),
189  };
190 
191  return ret;
192 }
193 
194 uint32_t ff_oklab_int_to_srgb_u8(struct Lab c)
195 {
196  const int64_t l_ = c.L + div_round64(25974LL * c.a, K) + div_round64(14143LL * c.b, K);
197  const int64_t m_ = c.L + div_round64(-6918LL * c.a, K) + div_round64(-4185LL * c.b, K);
198  const int64_t s_ = c.L + div_round64(-5864LL * c.a, K) + div_round64(-84638LL * c.b, K);
199 
200  const int32_t l = l_*l_*l_ / K2;
201  const int32_t m = m_*m_*m_ / K2;
202  const int32_t s = s_*s_*s_ / K2;
203 
204  const uint8_t r = ff_linear_int_to_srgb_u8((267169LL * l + -216771LL * m + 15137LL * s + K/2) / K);
205  const uint8_t g = ff_linear_int_to_srgb_u8((-83127LL * l + 171030LL * m + -22368LL * s + K/2) / K);
206  const uint8_t b = ff_linear_int_to_srgb_u8((-275LL * l + -46099LL * m + 111909LL * s + K/2) / K);
207 
208  return r<<16 | g<<8 | b;
209 }
210 
211 uint32_t ff_lowbias32(uint32_t x)
212 {
213  x ^= x >> 16;
214  x *= 0x7feb352d;
215  x ^= x >> 15;
216  x *= 0x846ca68b;
217  x ^= x >> 16;
218  return x;
219 }
r
const char * r
Definition: vf_curves.c:126
u
#define u(width, name, range_min, range_max)
Definition: cbs_h2645.c:262
palette.h
b
#define b
Definition: input.c:41
ff_oklab_int_to_srgb_u8
uint32_t ff_oklab_int_to_srgb_u8(struct Lab c)
OkLab to sRGB (non-linear) conversion.
Definition: palette.c:194
Lab::b
int32_t b
Definition: palette.h:31
ff_srgb_u8_to_linear_int
int32_t ff_srgb_u8_to_linear_int(uint8_t x)
Map sRGB 8-bit color component to a 16-bit linear value (gamma expand from electrical to optical valu...
Definition: palette.c:113
s
#define s(width, name)
Definition: cbs_vp9.c:256
g
const char * g
Definition: vf_curves.c:127
ff_lowbias32
uint32_t ff_lowbias32(uint32_t x)
Definition: palette.c:211
K2
#define K2
Definition: palette.c:26
c
Undefined Behavior In the C some operations are like signed integer dereferencing freed accessing outside allocated Undefined Behavior must not occur in a C it is not safe even if the output of undefined operations is unused The unsafety may seem nit picking but Optimizing compilers have in fact optimized code on the assumption that no undefined Behavior occurs Optimizing code based on wrong assumptions can and has in some cases lead to effects beyond the output of computations The signed integer overflow problem in speed critical code Code which is highly optimized and works with signed integers sometimes has the problem that often the output of the computation does not c
Definition: undefined.txt:32
Lab
Definition: palette.h:30
a
The reader does not expect b to be semantically here and if the code is changed by maybe adding a a division or other the signedness will almost certainly be mistaken To avoid this confusion a new type was SUINT is the C unsigned type but it holds a signed int to use the same example SUINT a
Definition: undefined.txt:41
i
#define i(width, name, range_min, range_max)
Definition: cbs_h2645.c:269
cbrt01_int
static int32_t cbrt01_int(int32_t x)
Definition: palette.c:135
common.h
linear2srgb
static const uint8_t linear2srgb[P+1]
Table mapping formula: f(x) = x < 0.0031308 ? x*12.92 : 1.055*x^(1/2.4)-0.055 (sRGB OETF) Where x is ...
Definition: palette.c:78
ret
ret
Definition: filter_design.txt:187
ff_linear_int_to_srgb_u8
uint8_t ff_linear_int_to_srgb_u8(int32_t x)
Map a 16-bit linear value to a sRGB 8-bit color component (gamma compressed from optical to electrica...
Definition: palette.c:118
K
#define K
Definition: palette.c:25
int32_t
int32_t
Definition: audioconvert.c:56
srgb2linear
static const uint16_t srgb2linear[256]
Table mapping formula: f(x) = x < 0.04045 ? x/12.92 : ((x+0.055)/1.055)^2.4 (sRGB EOTF) Where x is th...
Definition: palette.c:35
div_round64
static int64_t div_round64(int64_t a, int64_t b)
Definition: palette.c:168
P
#define P
Definition: palette.c:27
ff_srgb_u8_to_oklab_int
struct Lab ff_srgb_u8_to_oklab_int(uint32_t srgb)
sRGB (non-linear) to OkLab conversion
Definition: palette.c:170
Generated on Tue Feb 28 2023 21:33:51 for FFmpeg by   doxygen 1.8.17