/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 "./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)
static uint32_t kGammaToLinearTabS[GAMMA_TO_LINEAR_TAB_SIZE + 2];
#define LINEAR_TO_GAMMA_TAB_BITS 9
#define LINEAR_TO_GAMMA_TAB_SIZE (1 << LINEAR_TO_GAMMA_TAB_BITS)
static uint32_t kLinearToGammaTabS[LINEAR_TO_GAMMA_TAB_SIZE + 2];

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 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 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];
    }
    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, 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 ToLinearSrgb(uint16_t v, int bit_depth) {
  const int shift = GAMMA_TO_LINEAR_TAB_BITS - bit_depth;
  assert(v <= ((1 << bit_depth) - 1));
  if (shift > 0) {
    return kGammaToLinearTabS[v << shift];
  }
  return FixedPointInterpolation(v, kGammaToLinearTabS, -shift, 0);
}

static uint16_t FromLinearSrgb(uint32_t value, int bit_depth) {
  assert(value <= (1 << GAMMA_TO_LINEAR_BITS));
  return FixedPointInterpolation(
      value, kLinearToGammaTabS,
      (GAMMA_TO_LINEAR_BITS - LINEAR_TO_GAMMA_TAB_BITS),
      bit_depth - GAMMA_TO_LINEAR_BITS);
}

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

#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 ToLinear709(float gamma) {
  if (gamma < 0.f) {
    return 0.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 FromLinear709(float linear) {
  if (linear < 0.f) {
    return 0.f;
  } 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 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 == kSharpYuvTransferFunctionSrgb) {
    return ToLinearSrgb(v, bit_depth);
  }
  v_float = (float)v / ((1 << bit_depth) - 1);
  switch (transfer_type) {
    case kSharpYuvTransferFunctionBt709:
    case kSharpYuvTransferFunctionBt601:
    case kSharpYuvTransferFunctionBt2020_10Bit:
    case kSharpYuvTransferFunctionBt2020_12Bit:
      linear = ToLinear709(v_float);
      break;
    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 == kSharpYuvTransferFunctionSrgb) {
    return FromLinearSrgb(v, bit_depth);
  }
  v_float = (float)v / ((1 << 16) - 1);
  switch (transfer_type) {
    case kSharpYuvTransferFunctionBt709:
    case kSharpYuvTransferFunctionBt601:
    case kSharpYuvTransferFunctionBt2020_10Bit:
    case kSharpYuvTransferFunctionBt2020_12Bit:
      linear = FromLinear709(v_float);
      break;
    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-05-23 06:56

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