tor-browser

The Tor Browser
git clone https://git.dasho.dev/tor-browser.git
Log | Files | Refs | README | LICENSE

sharpyuv_gamma.c (13621B)

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