1 // Copyright 2022 Google Inc. All Rights Reserved.
2 //
3 // Use of this source code is governed by a BSD-style license
4 // that can be found in the COPYING file in the root of the source
5 // tree. An additional intellectual property rights grant can be found
6 // in the file PATENTS. All contributing project authors may
7 // be found in the AUTHORS file in the root of the source tree.
8 // -----------------------------------------------------------------------------
9 //
10 // Gamma correction utilities.
11
12 #include "sharpyuv/sharpyuv_gamma.h"
13
14 #include <assert.h>
15 #include <float.h>
16 #include <math.h>
17
18 #include "src/webp/types.h"
19
20 // Gamma correction compensates loss of resolution during chroma subsampling.
21 // Size of pre-computed table for converting from gamma to linear.
22 #define GAMMA_TO_LINEAR_TAB_BITS 10
23 #define GAMMA_TO_LINEAR_TAB_SIZE (1 << GAMMA_TO_LINEAR_TAB_BITS)
24 static uint32_t kGammaToLinearTabS[GAMMA_TO_LINEAR_TAB_SIZE + 2];
25 #define LINEAR_TO_GAMMA_TAB_BITS 9
26 #define LINEAR_TO_GAMMA_TAB_SIZE (1 << LINEAR_TO_GAMMA_TAB_BITS)
27 static uint32_t kLinearToGammaTabS[LINEAR_TO_GAMMA_TAB_SIZE + 2];
28
29 static const double kGammaF = 1. / 0.45;
30 #define GAMMA_TO_LINEAR_BITS 16
31
32 static volatile int kGammaTablesSOk = 0;
SharpYuvInitGammaTables(void)33 void SharpYuvInitGammaTables(void) {
34 assert(GAMMA_TO_LINEAR_BITS <= 16);
35 if (!kGammaTablesSOk) {
36 int v;
37 const double a = 0.09929682680944;
38 const double thresh = 0.018053968510807;
39 const double final_scale = 1 << GAMMA_TO_LINEAR_BITS;
40 // Precompute gamma to linear table.
41 {
42 const double norm = 1. / GAMMA_TO_LINEAR_TAB_SIZE;
43 const double a_rec = 1. / (1. + a);
44 for (v = 0; v <= GAMMA_TO_LINEAR_TAB_SIZE; ++v) {
45 const double g = norm * v;
46 double value;
47 if (g <= thresh * 4.5) {
48 value = g / 4.5;
49 } else {
50 value = pow(a_rec * (g + a), kGammaF);
51 }
52 kGammaToLinearTabS[v] = (uint32_t)(value * final_scale + .5);
53 }
54 // to prevent small rounding errors to cause read-overflow:
55 kGammaToLinearTabS[GAMMA_TO_LINEAR_TAB_SIZE + 1] =
56 kGammaToLinearTabS[GAMMA_TO_LINEAR_TAB_SIZE];
57 }
58 // Precompute linear to gamma table.
59 {
60 const double scale = 1. / LINEAR_TO_GAMMA_TAB_SIZE;
61 for (v = 0; v <= LINEAR_TO_GAMMA_TAB_SIZE; ++v) {
62 const double g = scale * v;
63 double value;
64 if (g <= thresh) {
65 value = 4.5 * g;
66 } else {
67 value = (1. + a) * pow(g, 1. / kGammaF) - a;
68 }
69 kLinearToGammaTabS[v] =
70 (uint32_t)(final_scale * value + 0.5);
71 }
72 // to prevent small rounding errors to cause read-overflow:
73 kLinearToGammaTabS[LINEAR_TO_GAMMA_TAB_SIZE + 1] =
74 kLinearToGammaTabS[LINEAR_TO_GAMMA_TAB_SIZE];
75 }
76 kGammaTablesSOk = 1;
77 }
78 }
79
Shift(int v,int shift)80 static WEBP_INLINE int Shift(int v, int shift) {
81 return (shift >= 0) ? (v << shift) : (v >> -shift);
82 }
83
FixedPointInterpolation(int v,uint32_t * tab,int tab_pos_shift_right,int tab_value_shift)84 static WEBP_INLINE uint32_t FixedPointInterpolation(int v, uint32_t* tab,
85 int tab_pos_shift_right,
86 int tab_value_shift) {
87 const uint32_t tab_pos = Shift(v, -tab_pos_shift_right);
88 // fractional part, in 'tab_pos_shift' fixed-point precision
89 const uint32_t x = v - (tab_pos << tab_pos_shift_right); // fractional part
90 // v0 / v1 are in kGammaToLinearBits fixed-point precision (range [0..1])
91 const uint32_t v0 = Shift(tab[tab_pos + 0], tab_value_shift);
92 const uint32_t v1 = Shift(tab[tab_pos + 1], tab_value_shift);
93 // Final interpolation.
94 const uint32_t v2 = (v1 - v0) * x; // note: v1 >= v0.
95 const int half =
96 (tab_pos_shift_right > 0) ? 1 << (tab_pos_shift_right - 1) : 0;
97 const uint32_t result = v0 + ((v2 + half) >> tab_pos_shift_right);
98 return result;
99 }
100
ToLinearSrgb(uint16_t v,int bit_depth)101 static uint32_t ToLinearSrgb(uint16_t v, int bit_depth) {
102 const int shift = GAMMA_TO_LINEAR_TAB_BITS - bit_depth;
103 if (shift > 0) {
104 return kGammaToLinearTabS[v << shift];
105 }
106 return FixedPointInterpolation(v, kGammaToLinearTabS, -shift, 0);
107 }
108
FromLinearSrgb(uint32_t value,int bit_depth)109 static uint16_t FromLinearSrgb(uint32_t value, int bit_depth) {
110 return FixedPointInterpolation(
111 value, kLinearToGammaTabS,
112 (GAMMA_TO_LINEAR_BITS - LINEAR_TO_GAMMA_TAB_BITS),
113 bit_depth - GAMMA_TO_LINEAR_BITS);
114 }
115
116 ////////////////////////////////////////////////////////////////////////////////
117
118 #define CLAMP(x, low, high) \
119 (((x) < (low)) ? (low) : (((high) < (x)) ? (high) : (x)))
120 #define MIN(a, b) (((a) < (b)) ? (a) : (b))
121 #define MAX(a, b) (((a) > (b)) ? (a) : (b))
122
Roundf(float x)123 static WEBP_INLINE float Roundf(float x) {
124 if (x < 0)
125 return (float)ceil((double)(x - 0.5f));
126 else
127 return (float)floor((double)(x + 0.5f));
128 }
129
Powf(float base,float exp)130 static WEBP_INLINE float Powf(float base, float exp) {
131 return (float)pow((double)base, (double)exp);
132 }
133
Log10f(float x)134 static WEBP_INLINE float Log10f(float x) { return (float)log10((double)x); }
135
ToLinear709(float gamma)136 static float ToLinear709(float gamma) {
137 if (gamma < 0.f) {
138 return 0.f;
139 } else if (gamma < 4.5f * 0.018053968510807f) {
140 return gamma / 4.5f;
141 } else if (gamma < 1.f) {
142 return Powf((gamma + 0.09929682680944f) / 1.09929682680944f, 1.f / 0.45f);
143 }
144 return 1.f;
145 }
146
FromLinear709(float linear)147 static float FromLinear709(float linear) {
148 if (linear < 0.f) {
149 return 0.f;
150 } else if (linear < 0.018053968510807f) {
151 return linear * 4.5f;
152 } else if (linear < 1.f) {
153 return 1.09929682680944f * Powf(linear, 0.45f) - 0.09929682680944f;
154 }
155 return 1.f;
156 }
157
ToLinear470M(float gamma)158 static float ToLinear470M(float gamma) {
159 return Powf(CLAMP(gamma, 0.f, 1.f), 2.2f);
160 }
161
FromLinear470M(float linear)162 static float FromLinear470M(float linear) {
163 return Powf(CLAMP(linear, 0.f, 1.f), 1.f / 2.2f);
164 }
165
ToLinear470Bg(float gamma)166 static float ToLinear470Bg(float gamma) {
167 return Powf(CLAMP(gamma, 0.f, 1.f), 2.8f);
168 }
169
FromLinear470Bg(float linear)170 static float FromLinear470Bg(float linear) {
171 return Powf(CLAMP(linear, 0.f, 1.f), 1.f / 2.8f);
172 }
173
ToLinearSmpte240(float gamma)174 static float ToLinearSmpte240(float gamma) {
175 if (gamma < 0.f) {
176 return 0.f;
177 } else if (gamma < 4.f * 0.022821585529445f) {
178 return gamma / 4.f;
179 } else if (gamma < 1.f) {
180 return Powf((gamma + 0.111572195921731f) / 1.111572195921731f, 1.f / 0.45f);
181 }
182 return 1.f;
183 }
184
FromLinearSmpte240(float linear)185 static float FromLinearSmpte240(float linear) {
186 if (linear < 0.f) {
187 return 0.f;
188 } else if (linear < 0.022821585529445f) {
189 return linear * 4.f;
190 } else if (linear < 1.f) {
191 return 1.111572195921731f * Powf(linear, 0.45f) - 0.111572195921731f;
192 }
193 return 1.f;
194 }
195
ToLinearLog100(float gamma)196 static float ToLinearLog100(float gamma) {
197 // The function is non-bijective so choose the middle of [0, 0.01].
198 const float mid_interval = 0.01f / 2.f;
199 return (gamma <= 0.0f) ? mid_interval
200 : Powf(10.0f, 2.f * (MIN(gamma, 1.f) - 1.0f));
201 }
202
FromLinearLog100(float linear)203 static float FromLinearLog100(float linear) {
204 return (linear < 0.01f) ? 0.0f : 1.0f + Log10f(MIN(linear, 1.f)) / 2.0f;
205 }
206
ToLinearLog100Sqrt10(float gamma)207 static float ToLinearLog100Sqrt10(float gamma) {
208 // The function is non-bijective so choose the middle of [0, 0.00316227766f[.
209 const float mid_interval = 0.00316227766f / 2.f;
210 return (gamma <= 0.0f) ? mid_interval
211 : Powf(10.0f, 2.5f * (MIN(gamma, 1.f) - 1.0f));
212 }
213
FromLinearLog100Sqrt10(float linear)214 static float FromLinearLog100Sqrt10(float linear) {
215 return (linear < 0.00316227766f) ? 0.0f
216 : 1.0f + Log10f(MIN(linear, 1.f)) / 2.5f;
217 }
218
ToLinearIec61966(float gamma)219 static float ToLinearIec61966(float gamma) {
220 if (gamma <= -4.5f * 0.018053968510807f) {
221 return Powf((-gamma + 0.09929682680944f) / -1.09929682680944f, 1.f / 0.45f);
222 } else if (gamma < 4.5f * 0.018053968510807f) {
223 return gamma / 4.5f;
224 }
225 return Powf((gamma + 0.09929682680944f) / 1.09929682680944f, 1.f / 0.45f);
226 }
227
FromLinearIec61966(float linear)228 static float FromLinearIec61966(float linear) {
229 if (linear <= -0.018053968510807f) {
230 return -1.09929682680944f * Powf(-linear, 0.45f) + 0.09929682680944f;
231 } else if (linear < 0.018053968510807f) {
232 return linear * 4.5f;
233 }
234 return 1.09929682680944f * Powf(linear, 0.45f) - 0.09929682680944f;
235 }
236
ToLinearBt1361(float gamma)237 static float ToLinearBt1361(float gamma) {
238 if (gamma < -0.25f) {
239 return -0.25f;
240 } else if (gamma < 0.f) {
241 return Powf((gamma - 0.02482420670236f) / -0.27482420670236f, 1.f / 0.45f) /
242 -4.f;
243 } else if (gamma < 4.5f * 0.018053968510807f) {
244 return gamma / 4.5f;
245 } else if (gamma < 1.f) {
246 return Powf((gamma + 0.09929682680944f) / 1.09929682680944f, 1.f / 0.45f);
247 }
248 return 1.f;
249 }
250
FromLinearBt1361(float linear)251 static float FromLinearBt1361(float linear) {
252 if (linear < -0.25f) {
253 return -0.25f;
254 } else if (linear < 0.f) {
255 return -0.27482420670236f * Powf(-4.f * linear, 0.45f) + 0.02482420670236f;
256 } else if (linear < 0.018053968510807f) {
257 return linear * 4.5f;
258 } else if (linear < 1.f) {
259 return 1.09929682680944f * Powf(linear, 0.45f) - 0.09929682680944f;
260 }
261 return 1.f;
262 }
263
ToLinearPq(float gamma)264 static float ToLinearPq(float gamma) {
265 if (gamma > 0.f) {
266 const float pow_gamma = Powf(gamma, 32.f / 2523.f);
267 const float num = MAX(pow_gamma - 107.f / 128.f, 0.0f);
268 const float den = MAX(2413.f / 128.f - 2392.f / 128.f * pow_gamma, FLT_MIN);
269 return Powf(num / den, 4096.f / 653.f);
270 }
271 return 0.f;
272 }
273
FromLinearPq(float linear)274 static float FromLinearPq(float linear) {
275 if (linear > 0.f) {
276 const float pow_linear = Powf(linear, 653.f / 4096.f);
277 const float num = 107.f / 128.f + 2413.f / 128.f * pow_linear;
278 const float den = 1.0f + 2392.f / 128.f * pow_linear;
279 return Powf(num / den, 2523.f / 32.f);
280 }
281 return 0.f;
282 }
283
ToLinearSmpte428(float gamma)284 static float ToLinearSmpte428(float gamma) {
285 return Powf(MAX(gamma, 0.f), 2.6f) / 0.91655527974030934f;
286 }
287
FromLinearSmpte428(float linear)288 static float FromLinearSmpte428(float linear) {
289 return Powf(0.91655527974030934f * MAX(linear, 0.f), 1.f / 2.6f);
290 }
291
292 // Conversion in BT.2100 requires RGB info. Simplify to gamma correction here.
ToLinearHlg(float gamma)293 static float ToLinearHlg(float gamma) {
294 if (gamma < 0.f) {
295 return 0.f;
296 } else if (gamma <= 0.5f) {
297 return Powf((gamma * gamma) * (1.f / 3.f), 1.2f);
298 }
299 return Powf((expf((gamma - 0.55991073f) / 0.17883277f) + 0.28466892f) / 12.0f,
300 1.2f);
301 }
302
FromLinearHlg(float linear)303 static float FromLinearHlg(float linear) {
304 linear = Powf(linear, 1.f / 1.2f);
305 if (linear < 0.f) {
306 return 0.f;
307 } else if (linear <= (1.f / 12.f)) {
308 return sqrtf(3.f * linear);
309 }
310 return 0.17883277f * logf(12.f * linear - 0.28466892f) + 0.55991073f;
311 }
312
SharpYuvGammaToLinear(uint16_t v,int bit_depth,SharpYuvTransferFunctionType transfer_type)313 uint32_t SharpYuvGammaToLinear(uint16_t v, int bit_depth,
314 SharpYuvTransferFunctionType transfer_type) {
315 float v_float, linear;
316 if (transfer_type == kSharpYuvTransferFunctionSrgb) {
317 return ToLinearSrgb(v, bit_depth);
318 }
319 v_float = (float)v / ((1 << bit_depth) - 1);
320 switch (transfer_type) {
321 case kSharpYuvTransferFunctionBt709:
322 case kSharpYuvTransferFunctionBt601:
323 case kSharpYuvTransferFunctionBt2020_10Bit:
324 case kSharpYuvTransferFunctionBt2020_12Bit:
325 linear = ToLinear709(v_float);
326 break;
327 case kSharpYuvTransferFunctionBt470M:
328 linear = ToLinear470M(v_float);
329 break;
330 case kSharpYuvTransferFunctionBt470Bg:
331 linear = ToLinear470Bg(v_float);
332 break;
333 case kSharpYuvTransferFunctionSmpte240:
334 linear = ToLinearSmpte240(v_float);
335 break;
336 case kSharpYuvTransferFunctionLinear:
337 return v;
338 case kSharpYuvTransferFunctionLog100:
339 linear = ToLinearLog100(v_float);
340 break;
341 case kSharpYuvTransferFunctionLog100_Sqrt10:
342 linear = ToLinearLog100Sqrt10(v_float);
343 break;
344 case kSharpYuvTransferFunctionIec61966:
345 linear = ToLinearIec61966(v_float);
346 break;
347 case kSharpYuvTransferFunctionBt1361:
348 linear = ToLinearBt1361(v_float);
349 break;
350 case kSharpYuvTransferFunctionSmpte2084:
351 linear = ToLinearPq(v_float);
352 break;
353 case kSharpYuvTransferFunctionSmpte428:
354 linear = ToLinearSmpte428(v_float);
355 break;
356 case kSharpYuvTransferFunctionHlg:
357 linear = ToLinearHlg(v_float);
358 break;
359 default:
360 assert(0);
361 linear = 0;
362 break;
363 }
364 return (uint32_t)Roundf(linear * ((1 << 16) - 1));
365 }
366
SharpYuvLinearToGamma(uint32_t v,int bit_depth,SharpYuvTransferFunctionType transfer_type)367 uint16_t SharpYuvLinearToGamma(uint32_t v, int bit_depth,
368 SharpYuvTransferFunctionType transfer_type) {
369 float v_float, linear;
370 if (transfer_type == kSharpYuvTransferFunctionSrgb) {
371 return FromLinearSrgb(v, bit_depth);
372 }
373 v_float = (float)v / ((1 << 16) - 1);
374 switch (transfer_type) {
375 case kSharpYuvTransferFunctionBt709:
376 case kSharpYuvTransferFunctionBt601:
377 case kSharpYuvTransferFunctionBt2020_10Bit:
378 case kSharpYuvTransferFunctionBt2020_12Bit:
379 linear = FromLinear709(v_float);
380 break;
381 case kSharpYuvTransferFunctionBt470M:
382 linear = FromLinear470M(v_float);
383 break;
384 case kSharpYuvTransferFunctionBt470Bg:
385 linear = FromLinear470Bg(v_float);
386 break;
387 case kSharpYuvTransferFunctionSmpte240:
388 linear = FromLinearSmpte240(v_float);
389 break;
390 case kSharpYuvTransferFunctionLinear:
391 return v;
392 case kSharpYuvTransferFunctionLog100:
393 linear = FromLinearLog100(v_float);
394 break;
395 case kSharpYuvTransferFunctionLog100_Sqrt10:
396 linear = FromLinearLog100Sqrt10(v_float);
397 break;
398 case kSharpYuvTransferFunctionIec61966:
399 linear = FromLinearIec61966(v_float);
400 break;
401 case kSharpYuvTransferFunctionBt1361:
402 linear = FromLinearBt1361(v_float);
403 break;
404 case kSharpYuvTransferFunctionSmpte2084:
405 linear = FromLinearPq(v_float);
406 break;
407 case kSharpYuvTransferFunctionSmpte428:
408 linear = FromLinearSmpte428(v_float);
409 break;
410 case kSharpYuvTransferFunctionHlg:
411 linear = FromLinearHlg(v_float);
412 break;
413 default:
414 assert(0);
415 linear = 0;
416 break;
417 }
418 return (uint16_t)Roundf(linear * ((1 << bit_depth) - 1));
419 }
420