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