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