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)
24static 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)
27static uint32_t kLinearToGammaTabS[LINEAR_TO_GAMMA_TAB_SIZE + 2];
28
29static const double kGammaF = 1. / 0.45;
30#define GAMMA_TO_LINEAR_BITS 16
31
32static volatile int kGammaTablesSOk = 0;
33void 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(x: a_rec * (g + a), y: 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(x: g, y: 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
80static WEBP_INLINE int Shift(int v, int shift) {
81 return (shift >= 0) ? (v << shift) : (v >> -shift);
82}
83
84static 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, shift: -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(v: tab[tab_pos + 0], shift: tab_value_shift);
92 const uint32_t v1 = Shift(v: tab[tab_pos + 1], shift: 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
101static 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, tab: kGammaToLinearTabS, tab_pos_shift_right: -shift, tab_value_shift: 0);
107}
108
109static uint16_t FromLinearSrgb(uint32_t value, int bit_depth) {
110 return FixedPointInterpolation(
111 v: value, tab: kLinearToGammaTabS,
112 tab_pos_shift_right: (GAMMA_TO_LINEAR_BITS - LINEAR_TO_GAMMA_TAB_BITS),
113 tab_value_shift: 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
123static WEBP_INLINE float Roundf(float x) {
124 if (x < 0)
125 return (float)ceil(x: (double)(x - 0.5f));
126 else
127 return (float)floor(x: (double)(x + 0.5f));
128}
129
130static WEBP_INLINE float Powf(float base, float exp) {
131 return (float)pow(x: (double)base, y: (double)exp);
132}
133
134static WEBP_INLINE float Log10f(float x) { return (float)log10(x: (double)x); }
135
136static 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(base: (gamma + 0.09929682680944f) / 1.09929682680944f, exp: 1.f / 0.45f);
143 }
144 return 1.f;
145}
146
147static 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(base: linear, exp: 0.45f) - 0.09929682680944f;
154 }
155 return 1.f;
156}
157
158static float ToLinear470M(float gamma) {
159 return Powf(CLAMP(gamma, 0.f, 1.f), exp: 2.2f);
160}
161
162static float FromLinear470M(float linear) {
163 return Powf(CLAMP(linear, 0.f, 1.f), exp: 1.f / 2.2f);
164}
165
166static float ToLinear470Bg(float gamma) {
167 return Powf(CLAMP(gamma, 0.f, 1.f), exp: 2.8f);
168}
169
170static float FromLinear470Bg(float linear) {
171 return Powf(CLAMP(linear, 0.f, 1.f), exp: 1.f / 2.8f);
172}
173
174static 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(base: (gamma + 0.111572195921731f) / 1.111572195921731f, exp: 1.f / 0.45f);
181 }
182 return 1.f;
183}
184
185static 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(base: linear, exp: 0.45f) - 0.111572195921731f;
192 }
193 return 1.f;
194}
195
196static 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(base: 10.0f, exp: 2.f * (MIN(gamma, 1.f) - 1.0f));
201}
202
203static float FromLinearLog100(float linear) {
204 return (linear < 0.01f) ? 0.0f : 1.0f + Log10f(MIN(linear, 1.f)) / 2.0f;
205}
206
207static 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(base: 10.0f, exp: 2.5f * (MIN(gamma, 1.f) - 1.0f));
212}
213
214static float FromLinearLog100Sqrt10(float linear) {
215 return (linear < 0.00316227766f) ? 0.0f
216 : 1.0f + Log10f(MIN(linear, 1.f)) / 2.5f;
217}
218
219static float ToLinearIec61966(float gamma) {
220 if (gamma <= -4.5f * 0.018053968510807f) {
221 return Powf(base: (-gamma + 0.09929682680944f) / -1.09929682680944f, exp: 1.f / 0.45f);
222 } else if (gamma < 4.5f * 0.018053968510807f) {
223 return gamma / 4.5f;
224 }
225 return Powf(base: (gamma + 0.09929682680944f) / 1.09929682680944f, exp: 1.f / 0.45f);
226}
227
228static float FromLinearIec61966(float linear) {
229 if (linear <= -0.018053968510807f) {
230 return -1.09929682680944f * Powf(base: -linear, exp: 0.45f) + 0.09929682680944f;
231 } else if (linear < 0.018053968510807f) {
232 return linear * 4.5f;
233 }
234 return 1.09929682680944f * Powf(base: linear, exp: 0.45f) - 0.09929682680944f;
235}
236
237static float ToLinearBt1361(float gamma) {
238 if (gamma < -0.25f) {
239 return -0.25f;
240 } else if (gamma < 0.f) {
241 return Powf(base: (gamma - 0.02482420670236f) / -0.27482420670236f, exp: 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(base: (gamma + 0.09929682680944f) / 1.09929682680944f, exp: 1.f / 0.45f);
247 }
248 return 1.f;
249}
250
251static float FromLinearBt1361(float linear) {
252 if (linear < -0.25f) {
253 return -0.25f;
254 } else if (linear < 0.f) {
255 return -0.27482420670236f * Powf(base: -4.f * linear, exp: 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(base: linear, exp: 0.45f) - 0.09929682680944f;
260 }
261 return 1.f;
262}
263
264static float ToLinearPq(float gamma) {
265 if (gamma > 0.f) {
266 const float pow_gamma = Powf(base: gamma, exp: 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(base: num / den, exp: 4096.f / 653.f);
270 }
271 return 0.f;
272}
273
274static float FromLinearPq(float linear) {
275 if (linear > 0.f) {
276 const float pow_linear = Powf(base: linear, exp: 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(base: num / den, exp: 2523.f / 32.f);
280 }
281 return 0.f;
282}
283
284static float ToLinearSmpte428(float gamma) {
285 return Powf(MAX(gamma, 0.f), exp: 2.6f) / 0.91655527974030934f;
286}
287
288static float FromLinearSmpte428(float linear) {
289 return Powf(base: 0.91655527974030934f * MAX(linear, 0.f), exp: 1.f / 2.6f);
290}
291
292// Conversion in BT.2100 requires RGB info. Simplify to gamma correction here.
293static float ToLinearHlg(float gamma) {
294 if (gamma < 0.f) {
295 return 0.f;
296 } else if (gamma <= 0.5f) {
297 return Powf(base: (gamma * gamma) * (1.f / 3.f), exp: 1.2f);
298 }
299 return Powf(base: (expf(x: (gamma - 0.55991073f) / 0.17883277f) + 0.28466892f) / 12.0f,
300 exp: 1.2f);
301}
302
303static float FromLinearHlg(float linear) {
304 linear = Powf(base: linear, exp: 1.f / 1.2f);
305 if (linear < 0.f) {
306 return 0.f;
307 } else if (linear <= (1.f / 12.f)) {
308 return sqrtf(x: 3.f * linear);
309 }
310 return 0.17883277f * logf(x: 12.f * linear - 0.28466892f) + 0.55991073f;
311}
312
313uint32_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(gamma: v_float);
326 break;
327 case kSharpYuvTransferFunctionBt470M:
328 linear = ToLinear470M(gamma: v_float);
329 break;
330 case kSharpYuvTransferFunctionBt470Bg:
331 linear = ToLinear470Bg(gamma: v_float);
332 break;
333 case kSharpYuvTransferFunctionSmpte240:
334 linear = ToLinearSmpte240(gamma: v_float);
335 break;
336 case kSharpYuvTransferFunctionLinear:
337 return v;
338 case kSharpYuvTransferFunctionLog100:
339 linear = ToLinearLog100(gamma: v_float);
340 break;
341 case kSharpYuvTransferFunctionLog100_Sqrt10:
342 linear = ToLinearLog100Sqrt10(gamma: v_float);
343 break;
344 case kSharpYuvTransferFunctionIec61966:
345 linear = ToLinearIec61966(gamma: v_float);
346 break;
347 case kSharpYuvTransferFunctionBt1361:
348 linear = ToLinearBt1361(gamma: v_float);
349 break;
350 case kSharpYuvTransferFunctionSmpte2084:
351 linear = ToLinearPq(gamma: v_float);
352 break;
353 case kSharpYuvTransferFunctionSmpte428:
354 linear = ToLinearSmpte428(gamma: v_float);
355 break;
356 case kSharpYuvTransferFunctionHlg:
357 linear = ToLinearHlg(gamma: v_float);
358 break;
359 default:
360 assert(0);
361 linear = 0;
362 break;
363 }
364 return (uint32_t)Roundf(x: linear * ((1 << 16) - 1));
365}
366
367uint16_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(value: 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(linear: v_float);
380 break;
381 case kSharpYuvTransferFunctionBt470M:
382 linear = FromLinear470M(linear: v_float);
383 break;
384 case kSharpYuvTransferFunctionBt470Bg:
385 linear = FromLinear470Bg(linear: v_float);
386 break;
387 case kSharpYuvTransferFunctionSmpte240:
388 linear = FromLinearSmpte240(linear: v_float);
389 break;
390 case kSharpYuvTransferFunctionLinear:
391 return v;
392 case kSharpYuvTransferFunctionLog100:
393 linear = FromLinearLog100(linear: v_float);
394 break;
395 case kSharpYuvTransferFunctionLog100_Sqrt10:
396 linear = FromLinearLog100Sqrt10(linear: v_float);
397 break;
398 case kSharpYuvTransferFunctionIec61966:
399 linear = FromLinearIec61966(linear: v_float);
400 break;
401 case kSharpYuvTransferFunctionBt1361:
402 linear = FromLinearBt1361(linear: v_float);
403 break;
404 case kSharpYuvTransferFunctionSmpte2084:
405 linear = FromLinearPq(linear: v_float);
406 break;
407 case kSharpYuvTransferFunctionSmpte428:
408 linear = FromLinearSmpte428(linear: v_float);
409 break;
410 case kSharpYuvTransferFunctionHlg:
411 linear = FromLinearHlg(linear: v_float);
412 break;
413 default:
414 assert(0);
415 linear = 0;
416 break;
417 }
418 return (uint16_t)Roundf(x: linear * ((1 << bit_depth) - 1));
419}
420

source code of qtimageformats/src/3rdparty/libwebp/sharpyuv/sharpyuv_gamma.c