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; |
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(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 | |
80 | static WEBP_INLINE int Shift(int v, int shift) { |
81 | return (shift >= 0) ? (v << shift) : (v >> -shift); |
82 | } |
83 | |
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, 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 | |
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, tab: kGammaToLinearTabS, tab_pos_shift_right: -shift, tab_value_shift: 0); |
107 | } |
108 | |
109 | static 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 | |
123 | static 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 | |
130 | static WEBP_INLINE float Powf(float base, float exp) { |
131 | return (float)pow(x: (double)base, y: (double)exp); |
132 | } |
133 | |
134 | static WEBP_INLINE float Log10f(float x) { return (float)log10(x: (double)x); } |
135 | |
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(base: (gamma + 0.09929682680944f) / 1.09929682680944f, exp: 1.f / 0.45f); |
143 | } |
144 | return 1.f; |
145 | } |
146 | |
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(base: linear, exp: 0.45f) - 0.09929682680944f; |
154 | } |
155 | return 1.f; |
156 | } |
157 | |
158 | static float ToLinear470M(float gamma) { |
159 | return Powf(CLAMP(gamma, 0.f, 1.f), exp: 2.2f); |
160 | } |
161 | |
162 | static float FromLinear470M(float linear) { |
163 | return Powf(CLAMP(linear, 0.f, 1.f), exp: 1.f / 2.2f); |
164 | } |
165 | |
166 | static float ToLinear470Bg(float gamma) { |
167 | return Powf(CLAMP(gamma, 0.f, 1.f), exp: 2.8f); |
168 | } |
169 | |
170 | static float FromLinear470Bg(float linear) { |
171 | return Powf(CLAMP(linear, 0.f, 1.f), exp: 1.f / 2.8f); |
172 | } |
173 | |
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(base: (gamma + 0.111572195921731f) / 1.111572195921731f, exp: 1.f / 0.45f); |
181 | } |
182 | return 1.f; |
183 | } |
184 | |
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(base: linear, exp: 0.45f) - 0.111572195921731f; |
192 | } |
193 | return 1.f; |
194 | } |
195 | |
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(base: 10.0f, exp: 2.f * (MIN(gamma, 1.f) - 1.0f)); |
201 | } |
202 | |
203 | static float FromLinearLog100(float linear) { |
204 | return (linear < 0.01f) ? 0.0f : 1.0f + Log10f(MIN(linear, 1.f)) / 2.0f; |
205 | } |
206 | |
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(base: 10.0f, exp: 2.5f * (MIN(gamma, 1.f) - 1.0f)); |
212 | } |
213 | |
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 | |
219 | static 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 | |
228 | static 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 | |
237 | static 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 | |
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(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 | |
264 | static 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 | |
274 | static 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 | |
284 | static float ToLinearSmpte428(float gamma) { |
285 | return Powf(MAX(gamma, 0.f), exp: 2.6f) / 0.91655527974030934f; |
286 | } |
287 | |
288 | static 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. |
293 | static 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 | |
303 | static 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 | |
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(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 | |
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(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 |
Definitions
- kGammaToLinearTabS
- kLinearToGammaTabS
- kGammaF
- kGammaTablesSOk
- SharpYuvInitGammaTables
- Shift
- FixedPointInterpolation
- ToLinearSrgb
- FromLinearSrgb
- Roundf
- Powf
- Log10f
- ToLinear709
- FromLinear709
- ToLinear470M
- FromLinear470M
- ToLinear470Bg
- FromLinear470Bg
- ToLinearSmpte240
- FromLinearSmpte240
- ToLinearLog100
- FromLinearLog100
- ToLinearLog100Sqrt10
- FromLinearLog100Sqrt10
- ToLinearIec61966
- FromLinearIec61966
- ToLinearBt1361
- FromLinearBt1361
- ToLinearPq
- FromLinearPq
- ToLinearSmpte428
- FromLinearSmpte428
- ToLinearHlg
- FromLinearHlg
- SharpYuvGammaToLinear
Learn to use CMake with our Intro Training
Find out more