media/libwebp/sharpyuv/sharpyuv_gamma.c

   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 <math.h>
  16
  17 #include "src/webp/types.h"
  18
  19 // Gamma correction compensates loss of resolution during chroma subsampling.
  20 // Size of pre-computed table for converting from gamma to linear.
  21 #define GAMMA_TO_LINEAR_TAB_BITS 10
  22 #define GAMMA_TO_LINEAR_TAB_SIZE (1 << GAMMA_TO_LINEAR_TAB_BITS)
  23 static uint32_t kGammaToLinearTabS[GAMMA_TO_LINEAR_TAB_SIZE + 2];
  24 #define LINEAR_TO_GAMMA_TAB_BITS 9
  25 #define LINEAR_TO_GAMMA_TAB_SIZE (1 << LINEAR_TO_GAMMA_TAB_BITS)
  26 static uint32_t kLinearToGammaTabS[LINEAR_TO_GAMMA_TAB_SIZE + 2];
  27
  28 static const double kGammaF = 1. / 0.45;
  29 #define GAMMA_TO_LINEAR_BITS 16
  30
  31 static volatile int kGammaTablesSOk = 0;
  32 void SharpYuvInitGammaTables(void) {
  33   assert(GAMMA_TO_LINEAR_BITS <= 16);
  34   if (!kGammaTablesSOk) {
  35     int v;
  36     const double a = 0.09929682680944;
  37     const double thresh = 0.018053968510807;
  38     const double final_scale = 1 << GAMMA_TO_LINEAR_BITS;
  39     // Precompute gamma to linear table.
  40     {
  41       const double norm = 1. / GAMMA_TO_LINEAR_TAB_SIZE;
  42       const double a_rec = 1. / (1. + a);
  43       for (v = 0; v <= GAMMA_TO_LINEAR_TAB_SIZE; ++v) {
  44         const double g = norm * v;
  45         double value;
  46         if (g <= thresh * 4.5) {
  47           value = g / 4.5;
  48         } else {
  49           value = pow(a_rec * (g + a), kGammaF);
  50         }
  51         kGammaToLinearTabS[v] = (uint32_t)(value * final_scale + .5);
  52       }
  53       // to prevent small rounding errors to cause read-overflow:
  54       kGammaToLinearTabS[GAMMA_TO_LINEAR_TAB_SIZE + 1] =
  55           kGammaToLinearTabS[GAMMA_TO_LINEAR_TAB_SIZE];
  56     }
  57     // Precompute linear to gamma table.
  58     {
  59       const double scale = 1. / LINEAR_TO_GAMMA_TAB_SIZE;
  60       for (v = 0; v <= LINEAR_TO_GAMMA_TAB_SIZE; ++v) {
  61         const double g = scale * v;
  62         double value;
  63         if (g <= thresh) {
  64           value = 4.5 * g;
  65         } else {
  66           value = (1. + a) * pow(g, 1. / kGammaF) - a;
  67         }
  68         kLinearToGammaTabS[v] =
  69             (uint32_t)(final_scale * value + 0.5);
  70       }
  71       // to prevent small rounding errors to cause read-overflow:
  72       kLinearToGammaTabS[LINEAR_TO_GAMMA_TAB_SIZE + 1] =
  73           kLinearToGammaTabS[LINEAR_TO_GAMMA_TAB_SIZE];
  74     }
  75     kGammaTablesSOk = 1;
  76   }
  77 }
  78
  79 static WEBP_INLINE int Shift(int v, int shift) {
  80   return (shift >= 0) ? (v << shift) : (v >> -shift);
  81 }
  82
  83 static WEBP_INLINE uint32_t FixedPointInterpolation(int v, uint32_t* tab,
  84                                                     int tab_pos_shift_right,
  85                                                     int tab_value_shift) {
  86   const uint32_t tab_pos = Shift(v, -tab_pos_shift_right);
  87   // fractional part, in 'tab_pos_shift' fixed-point precision
  88   const uint32_t x = v - (tab_pos << tab_pos_shift_right);  // fractional part
  89   // v0 / v1 are in kGammaToLinearBits fixed-point precision (range [0..1])
  90   const uint32_t v0 = Shift(tab[tab_pos + 0], tab_value_shift);
  91   const uint32_t v1 = Shift(tab[tab_pos + 1], tab_value_shift);
  92   // Final interpolation.
  93   const uint32_t v2 = (v1 - v0) * x;  // note: v1 >= v0.
  94   const int half =
  95       (tab_pos_shift_right > 0) ? 1 << (tab_pos_shift_right - 1) : 0;
  96   const uint32_t result = v0 + ((v2 + half) >> tab_pos_shift_right);
  97   return result;
  98 }
  99
 100 uint32_t SharpYuvGammaToLinear(uint16_t v, int bit_depth) {
 101   const int shift = GAMMA_TO_LINEAR_TAB_BITS - bit_depth;
 102   if (shift > 0) {
 103     return kGammaToLinearTabS[v << shift];
 104   }
 105   return FixedPointInterpolation(v, kGammaToLinearTabS, -shift, 0);
 106 }
 107
 108 uint16_t SharpYuvLinearToGamma(uint32_t value, int bit_depth) {
 109   return FixedPointInterpolation(
 110       value, kLinearToGammaTabS,
 111       (GAMMA_TO_LINEAR_BITS - LINEAR_TO_GAMMA_TAB_BITS),
 112       bit_depth - GAMMA_TO_LINEAR_BITS);
 113 }