/src/libwebp/sharpyuv/sharpyuv_gamma.c

Source
// Copyright 2022 Google Inc. All Rights Reserved.
//
// Use of this source code is governed by a BSD-style license
// that can be found in the COPYING file in the root of the source
// tree. An additional intellectual property rights grant can be found
// in the file PATENTS. All contributing project authors may
// be found in the AUTHORS file in the root of the source tree.
// -----------------------------------------------------------------------------
//
// Gamma correction utilities.

#include "./sharpyuv_gamma.h"

#include <assert.h>
#include <float.h>
#include <math.h>
#include <stdint.h>

#include "./sharpyuv.h"
#include "webp/types.h"

// Gamma correction compensates loss of resolution during chroma subsampling.
// Size of pre-computed table for converting from gamma to linear.
#define GAMMA_TO_LINEAR_TAB_BITS 10
#define GAMMA_TO_LINEAR_TAB_SIZE (1 << GAMMA_TO_LINEAR_TAB_BITS)
uint32_t kSharpYuvGammaToLinearTabS[GAMMA_TO_LINEAR_TAB_SIZE + 2];
#define kGammaToLinearTabS kSharpYuvGammaToLinearTabS
uint32_t kSharpYuvGammaToLinearTabSrgb[GAMMA_TO_LINEAR_TAB_SIZE + 2];
#define kGammaToLinearTabSrgb kSharpYuvGammaToLinearTabSrgb
#define LINEAR_TO_GAMMA_TAB_BITS 9
#define LINEAR_TO_GAMMA_TAB_SIZE (1 << LINEAR_TO_GAMMA_TAB_BITS)
uint32_t kSharpYuvLinearToGammaTabS[LINEAR_TO_GAMMA_TAB_SIZE + 2];
#define kLinearToGammaTabS kSharpYuvLinearToGammaTabS
uint32_t kSharpYuvLinearToGammaTabSrgb[LINEAR_TO_GAMMA_TAB_SIZE + 2];
#define kLinearToGammaTabSrgb kSharpYuvLinearToGammaTabSrgb

static const double kGammaF = 1. / 0.45;
#define GAMMA_TO_LINEAR_BITS 16

static volatile int kGammaTablesSOk = 0;
void SharpYuvInitGammaTables(void) {
  assert(GAMMA_TO_LINEAR_BITS <= 16);
  if (!kGammaTablesSOk) {
    int v;
    const double a = 0.09929682680944;
    const double thresh = 0.018053968510807;
    const double final_scale = 1 << GAMMA_TO_LINEAR_BITS;
    // Precompute BT.709 gamma to linear table.
    {
      const double norm = 1. / GAMMA_TO_LINEAR_TAB_SIZE;
      const double a_rec = 1. / (1. + a);
      for (v = 0; v <= GAMMA_TO_LINEAR_TAB_SIZE; ++v) {
        const double g = norm * v;
        double value;
        if (g <= thresh * 4.5) {
          value = g / 4.5;
        } else {
          value = pow(a_rec * (g + a), kGammaF);
        }
        kGammaToLinearTabS[v] = (uint32_t)(value * final_scale + .5);
      }
      // to prevent small rounding errors to cause read-overflow:
      kGammaToLinearTabS[GAMMA_TO_LINEAR_TAB_SIZE + 1] =
          kGammaToLinearTabS[GAMMA_TO_LINEAR_TAB_SIZE];
    }
    // Precompute BT.709 linear to gamma table.
    {
      const double scale = 1. / LINEAR_TO_GAMMA_TAB_SIZE;
      for (v = 0; v <= LINEAR_TO_GAMMA_TAB_SIZE; ++v) {
        const double g = scale * v;
        double value;
        if (g <= thresh) {
          value = 4.5 * g;
        } else {
          value = (1. + a) * pow(g, 1. / kGammaF) - a;
        }
        kLinearToGammaTabS[v] = (uint32_t)(final_scale * value + 0.5);
      }
      // to prevent small rounding errors to cause read-overflow:
      kLinearToGammaTabS[LINEAR_TO_GAMMA_TAB_SIZE + 1] =
          kLinearToGammaTabS[LINEAR_TO_GAMMA_TAB_SIZE];
    }
    // Precompute sRGB gamma to linear table (IEC 61966-2-1).
    {
      const double norm = 1. / GAMMA_TO_LINEAR_TAB_SIZE;
      for (v = 0; v <= GAMMA_TO_LINEAR_TAB_SIZE; ++v) {
        const double g = norm * v;
        double value;
        if (g <= 0.04045) {
          value = g / 12.92;
        } else {
          value = pow((g + 0.055) / 1.055, 2.4);
        }
        kGammaToLinearTabSrgb[v] = (uint32_t)(value * final_scale + .5);
      }
      kGammaToLinearTabSrgb[GAMMA_TO_LINEAR_TAB_SIZE + 1] =
          kGammaToLinearTabSrgb[GAMMA_TO_LINEAR_TAB_SIZE];
    }
    // Precompute sRGB linear to gamma table (IEC 61966-2-1).
    {
      const double scale = 1. / LINEAR_TO_GAMMA_TAB_SIZE;
      for (v = 0; v <= LINEAR_TO_GAMMA_TAB_SIZE; ++v) {
        const double g = scale * v;
        double value;
        if (g <= 0.0031308) {
          value = 12.92 * g;
        } else {
          value = 1.055 * pow(g, 1. / 2.4) - 0.055;
        }
        kLinearToGammaTabSrgb[v] = (uint32_t)(final_scale * value + 0.5);
      }
      kLinearToGammaTabSrgb[LINEAR_TO_GAMMA_TAB_SIZE + 1] =
          kLinearToGammaTabSrgb[LINEAR_TO_GAMMA_TAB_SIZE];
    }
    kGammaTablesSOk = 1;
  }
}

static WEBP_INLINE int Shift(int v, int shift) {
  return (shift >= 0) ? (v << shift) : (v >> -shift);
}

static WEBP_INLINE uint32_t FixedPointInterpolation(int v, const uint32_t* tab,
                                                    int tab_pos_shift_right,
                                                    int tab_value_shift) {
  const uint32_t tab_pos = Shift(v, -tab_pos_shift_right);
  // fractional part, in 'tab_pos_shift' fixed-point precision
  const uint32_t x = v - (tab_pos << tab_pos_shift_right);  // fractional part
  // v0 / v1 are in kGammaToLinearBits fixed-point precision (range [0..1])
  const uint32_t v0 = Shift(tab[tab_pos + 0], tab_value_shift);
  const uint32_t v1 = Shift(tab[tab_pos + 1], tab_value_shift);
  // Final interpolation.
  const uint32_t v2 = (v1 - v0) * x;  // note: v1 >= v0.
  const int half =
      (tab_pos_shift_right > 0) ? 1 << (tab_pos_shift_right - 1) : 0;
  const uint32_t result = v0 + ((v2 + half) >> tab_pos_shift_right);
  return result;
}

static uint32_t ToLinear(uint16_t v, int bit_depth, const uint32_t* tab) {
  const int shift = GAMMA_TO_LINEAR_TAB_BITS - bit_depth;
  assert(v <= ((1 << bit_depth) - 1));
  if (shift >= 0) {
    // shift == 0 is a direct lookup, no need for interpolation
    return tab[v << shift];
  }
  return FixedPointInterpolation(v, tab, -shift, 0);
}

static uint16_t FromLinear(uint32_t value, int bit_depth, const uint32_t* tab) {
  const uint32_t max = (1 << bit_depth) - 1;
  uint32_t v;
  assert(value <= (1 << GAMMA_TO_LINEAR_BITS));
  v = FixedPointInterpolation(value, tab,
                              (GAMMA_TO_LINEAR_BITS - LINEAR_TO_GAMMA_TAB_BITS),
                              bit_depth - GAMMA_TO_LINEAR_BITS);
  return (v > max) ? (uint16_t)max : (uint16_t)v;
}

////////////////////////////////////////////////////////////////////////////////

#define CLAMP(x, low, high) \
  (((x) < (low)) ? (low) : (((high) < (x)) ? (high) : (x)))
#define MIN(a, b) (((a) < (b)) ? (a) : (b))
#define MAX(a, b) (((a) > (b)) ? (a) : (b))

static WEBP_INLINE float Roundf(float x) {
  if (x < 0) {
    return (float)ceil((double)(x - 0.5f));
  } else {
    return (float)floor((double)(x + 0.5f));
  }
}

static WEBP_INLINE float Powf(float base, float exp) {
  return (float)pow((double)base, (double)exp);
}

static WEBP_INLINE float Log10f(float x) { return (float)log10((double)x); }

static float ToLinear470M(float gamma) {
  return Powf(CLAMP(gamma, 0.f, 1.f), 2.2f);
}

static float FromLinear470M(float linear) {
  return Powf(CLAMP(linear, 0.f, 1.f), 1.f / 2.2f);
}

static float ToLinear470Bg(float gamma) {
  return Powf(CLAMP(gamma, 0.f, 1.f), 2.8f);
}

static float FromLinear470Bg(float linear) {
  return Powf(CLAMP(linear, 0.f, 1.f), 1.f / 2.8f);
}

static float ToLinearSmpte240(float gamma) {
  if (gamma < 0.f) {
    return 0.f;
  } else if (gamma < 4.f * 0.022821585529445f) {
    return gamma / 4.f;
  } else if (gamma < 1.f) {
    return Powf((gamma + 0.111572195921731f) / 1.111572195921731f, 1.f / 0.45f);
  }
  return 1.f;
}

static float FromLinearSmpte240(float linear) {
  if (linear < 0.f) {
    return 0.f;
  } else if (linear < 0.022821585529445f) {
    return linear * 4.f;
  } else if (linear < 1.f) {
    return 1.111572195921731f * Powf(linear, 0.45f) - 0.111572195921731f;
  }
  return 1.f;
}

static float ToLinearLog100(float gamma) {
  // The function is non-bijective so choose the middle of [0, 0.01].
  const float mid_interval = 0.01f / 2.f;
  return (gamma <= 0.0f) ? mid_interval
                         : Powf(10.0f, 2.f * (MIN(gamma, 1.f) - 1.0f));
}

static float FromLinearLog100(float linear) {
  return (linear < 0.01f) ? 0.0f : 1.0f + Log10f(MIN(linear, 1.f)) / 2.0f;
}

static float ToLinearLog100Sqrt10(float gamma) {
  // The function is non-bijective so choose the middle of [0, 0.00316227766f[.
  const float mid_interval = 0.00316227766f / 2.f;
  return (gamma <= 0.0f) ? mid_interval
                         : Powf(10.0f, 2.5f * (MIN(gamma, 1.f) - 1.0f));
}

static float FromLinearLog100Sqrt10(float linear) {
  return (linear < 0.00316227766f) ? 0.0f
                                   : 1.0f + Log10f(MIN(linear, 1.f)) / 2.5f;
}

static float ToLinearIec61966(float gamma) {
  if (gamma <= -4.5f * 0.018053968510807f) {
    return Powf((-gamma + 0.09929682680944f) / -1.09929682680944f, 1.f / 0.45f);
  } else if (gamma < 4.5f * 0.018053968510807f) {
    return gamma / 4.5f;
  }
  return Powf((gamma + 0.09929682680944f) / 1.09929682680944f, 1.f / 0.45f);
}

static float FromLinearIec61966(float linear) {
  if (linear <= -0.018053968510807f) {
    return -1.09929682680944f * Powf(-linear, 0.45f) + 0.09929682680944f;
  } else if (linear < 0.018053968510807f) {
    return linear * 4.5f;
  }
  return 1.09929682680944f * Powf(linear, 0.45f) - 0.09929682680944f;
}

static float ToLinearBt1361(float gamma) {
  if (gamma < -0.25f) {
    return -0.25f;
  } else if (gamma < 0.f) {
    return Powf((gamma - 0.02482420670236f) / -0.27482420670236f, 1.f / 0.45f) /
           -4.f;
  } else if (gamma < 4.5f * 0.018053968510807f) {
    return gamma / 4.5f;
  } else if (gamma < 1.f) {
    return Powf((gamma + 0.09929682680944f) / 1.09929682680944f, 1.f / 0.45f);
  }
  return 1.f;
}

static float FromLinearBt1361(float linear) {
  if (linear < -0.25f) {
    return -0.25f;
  } else if (linear < 0.f) {
    return -0.27482420670236f * Powf(-4.f * linear, 0.45f) + 0.02482420670236f;
  } else if (linear < 0.018053968510807f) {
    return linear * 4.5f;
  } else if (linear < 1.f) {
    return 1.09929682680944f * Powf(linear, 0.45f) - 0.09929682680944f;
  }
  return 1.f;
}

static float ToLinearPq(float gamma) {
  if (gamma > 0.f) {
    const float pow_gamma = Powf(gamma, 32.f / 2523.f);
    const float num = MAX(pow_gamma - 107.f / 128.f, 0.0f);
    const float den = MAX(2413.f / 128.f - 2392.f / 128.f * pow_gamma, FLT_MIN);
    return Powf(num / den, 4096.f / 653.f);
  }
  return 0.f;
}

static float FromLinearPq(float linear) {
  if (linear > 0.f) {
    const float pow_linear = Powf(linear, 653.f / 4096.f);
    const float num = 107.f / 128.f + 2413.f / 128.f * pow_linear;
    const float den = 1.0f + 2392.f / 128.f * pow_linear;
    return Powf(num / den, 2523.f / 32.f);
  }
  return 0.f;
}

static float ToLinearSmpte428(float gamma) {
  return Powf(MAX(gamma, 0.f), 2.6f) / 0.91655527974030934f;
}

static float FromLinearSmpte428(float linear) {
  return Powf(0.91655527974030934f * MAX(linear, 0.f), 1.f / 2.6f);
}

// Conversion in BT.2100 requires RGB info. Simplify to gamma correction here.
static float ToLinearHlg(float gamma) {
  if (gamma < 0.f) {
    return 0.f;
  } else if (gamma <= 0.5f) {
    return Powf((gamma * gamma) * (1.f / 3.f), 1.2f);
  }
  return Powf((expf((gamma - 0.55991073f) / 0.17883277f) + 0.28466892f) / 12.0f,
              1.2f);
}

static float FromLinearHlg(float linear) {
  linear = Powf(linear, 1.f / 1.2f);
  if (linear < 0.f) {
    return 0.f;
  } else if (linear <= (1.f / 12.f)) {
    return sqrtf(3.f * linear);
  }
  return 0.17883277f * logf(12.f * linear - 0.28466892f) + 0.55991073f;
}

uint32_t SharpYuvGammaToLinear(uint16_t v, int bit_depth,
                               SharpYuvTransferFunctionType transfer_type) {
  float v_float, linear;
  if (transfer_type == kSharpYuvTransferFunctionBt709 ||
      transfer_type == kSharpYuvTransferFunctionBt601 ||
      transfer_type == kSharpYuvTransferFunctionBt2020_10Bit ||
      transfer_type == kSharpYuvTransferFunctionBt2020_12Bit) {
    return ToLinear(v, bit_depth, kGammaToLinearTabS);
  }
  if (transfer_type == kSharpYuvTransferFunctionSrgb) {
    return ToLinear(v, bit_depth, kGammaToLinearTabSrgb);
  }
  v_float = (float)v / ((1 << bit_depth) - 1);
  switch (transfer_type) {
    case kSharpYuvTransferFunctionBt470M:
      linear = ToLinear470M(v_float);
      break;
    case kSharpYuvTransferFunctionBt470Bg:
      linear = ToLinear470Bg(v_float);
      break;
    case kSharpYuvTransferFunctionSmpte240:
      linear = ToLinearSmpte240(v_float);
      break;
    case kSharpYuvTransferFunctionLinear:
      return v;
    case kSharpYuvTransferFunctionLog100:
      linear = ToLinearLog100(v_float);
      break;
    case kSharpYuvTransferFunctionLog100_Sqrt10:
      linear = ToLinearLog100Sqrt10(v_float);
      break;
    case kSharpYuvTransferFunctionIec61966:
      linear = ToLinearIec61966(v_float);
      break;
    case kSharpYuvTransferFunctionBt1361:
      linear = ToLinearBt1361(v_float);
      break;
    case kSharpYuvTransferFunctionSmpte2084:
      linear = ToLinearPq(v_float);
      break;
    case kSharpYuvTransferFunctionSmpte428:
      linear = ToLinearSmpte428(v_float);
      break;
    case kSharpYuvTransferFunctionHlg:
      linear = ToLinearHlg(v_float);
      break;
    default:
      assert(0);
      linear = 0;
      break;
  }
  return (uint32_t)Roundf(linear * ((1 << 16) - 1));
}

uint16_t SharpYuvLinearToGamma(uint32_t v, int bit_depth,
                               SharpYuvTransferFunctionType transfer_type) {
  float v_float, linear;
  if (transfer_type == kSharpYuvTransferFunctionBt709 ||
      transfer_type == kSharpYuvTransferFunctionBt601 ||
      transfer_type == kSharpYuvTransferFunctionBt2020_10Bit ||
      transfer_type == kSharpYuvTransferFunctionBt2020_12Bit) {
    return FromLinear(v, bit_depth, kLinearToGammaTabS);
  }
  if (transfer_type == kSharpYuvTransferFunctionSrgb) {
    return FromLinear(v, bit_depth, kLinearToGammaTabSrgb);
  }
  v_float = (float)v / ((1 << 16) - 1);
  switch (transfer_type) {
    case kSharpYuvTransferFunctionBt470M:
      linear = FromLinear470M(v_float);
      break;
    case kSharpYuvTransferFunctionBt470Bg:
      linear = FromLinear470Bg(v_float);
      break;
    case kSharpYuvTransferFunctionSmpte240:
      linear = FromLinearSmpte240(v_float);
      break;
    case kSharpYuvTransferFunctionLinear:
      return v;
    case kSharpYuvTransferFunctionLog100:
      linear = FromLinearLog100(v_float);
      break;
    case kSharpYuvTransferFunctionLog100_Sqrt10:
      linear = FromLinearLog100Sqrt10(v_float);
      break;
    case kSharpYuvTransferFunctionIec61966:
      linear = FromLinearIec61966(v_float);
      break;
    case kSharpYuvTransferFunctionBt1361:
      linear = FromLinearBt1361(v_float);
      break;
    case kSharpYuvTransferFunctionSmpte2084:
      linear = FromLinearPq(v_float);
      break;
    case kSharpYuvTransferFunctionSmpte428:
      linear = FromLinearSmpte428(v_float);
      break;
    case kSharpYuvTransferFunctionHlg:
      linear = FromLinearHlg(v_float);
      break;
    default:
      assert(0);
      linear = 0;
      break;
  }
  return (uint16_t)Roundf(linear * ((1 << bit_depth) - 1));
}

Coverage Report

Created: 2026-11-12 07:35

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