Coverage Report

Created: 2026-09-28 07:02

next uncovered line (L), next uncovered region (R), next uncovered branch (B)
/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
}