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 }