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