1 /* 2 * This file is part of RawTherapee. 3 * 4 * Copyright (c) 2004-2010 Gabor Horvath <hgabor@rawtherapee.com> 5 * 6 * RawTherapee is free software: you can redistribute it and/or modify 7 * it under the terms of the GNU General Public License as published by 8 * the Free Software Foundation, either version 3 of the License, or 9 * (at your option) any later version. 10 * 11 * RawTherapee is distributed in the hope that it will be useful, 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 * GNU General Public License for more details. 15 * 16 * You should have received a copy of the GNU General Public License 17 * along with RawTherapee. If not, see <https://www.gnu.org/licenses/>. 18 */ 19 20 #pragma once 21 22 #include <array> 23 #include <glibmm/ustring.h> 24 25 #include "rt_math.h" 26 #include "LUT.h" 27 #include "iccmatrices.h" 28 #include "lcms2.h" 29 #include "sleef.h" 30 31 #define SAT(a,b,c) ((float)max(a,b,c)-(float)min(a,b,c))/(float)max(a,b,c) 32 33 namespace rtengine 34 { 35 36 typedef std::array<double, 7> GammaValues; 37 38 #ifdef _DEBUG 39 40 class MunsellDebugInfo 41 { 42 public: 43 float maxdhuelum[4]; 44 float maxdhue[4]; 45 unsigned int depass; 46 unsigned int depassLum; 47 48 MunsellDebugInfo(); 49 void reinitValues(); 50 }; 51 52 #endif 53 54 55 class Color 56 { 57 58 private: 59 // Jacques' 195 LUTf for Munsell Lch correction 60 static LUTf _4P10, _4P20, _4P30, _4P40, _4P50, _4P60; 61 static LUTf _1P10, _1P20, _1P30, _1P40, _1P50, _1P60; 62 static LUTf _5B40, _5B50, _5B60, _5B70, _5B80; 63 static LUTf _7B40, _7B50, _7B60, _7B70, _7B80; 64 static LUTf _9B40, _9B50, _9B60, _9B70, _9B80; 65 static LUTf _10B40, _10B50, _10B60, _10B70, _10B80; 66 static LUTf _05PB40, _05PB50, _05PB60, _05PB70, _05PB80; 67 static LUTf _10PB10, _10PB20, _10PB30, _10PB40, _10PB50, _10PB60; 68 static LUTf _9PB10, _9PB20, _9PB30, _9PB40, _9PB50, _9PB60, _9PB70, _9PB80; 69 static LUTf _75PB10, _75PB20, _75PB30, _75PB40, _75PB50, _75PB60, _75PB70, _75PB80; 70 static LUTf _6PB10, _6PB20, _6PB30, _6PB40, _6PB50, _6PB60, _6PB70, _6PB80; 71 static LUTf _45PB10, _45PB20, _45PB30, _45PB40, _45PB50, _45PB60, _45PB70, _45PB80; 72 static LUTf _3PB10, _3PB20, _3PB30, _3PB40, _3PB50, _3PB60, _3PB70, _3PB80; 73 static LUTf _15PB10, _15PB20, _15PB30, _15PB40, _15PB50, _15PB60, _15PB70, _15PB80; 74 static LUTf _10YR20, _10YR30, _10YR40, _10YR50, _10YR60, _10YR70, _10YR80, _10YR90; 75 static LUTf _85YR20, _85YR30, _85YR40, _85YR50, _85YR60, _85YR70, _85YR80, _85YR90; 76 static LUTf _7YR30, _7YR40, _7YR50, _7YR60, _7YR70, _7YR80; 77 static LUTf _55YR30, _55YR40, _55YR50, _55YR60, _55YR70, _55YR80, _55YR90; 78 static LUTf _4YR30, _4YR40, _4YR50, _4YR60, _4YR70, _4YR80; 79 static LUTf _25YR30, _25YR40, _25YR50, _25YR60, _25YR70; 80 static LUTf _10R30, _10R40, _10R50, _10R60, _10R70; 81 static LUTf _9R30, _9R40, _9R50, _9R60, _9R70; 82 static LUTf _7R30, _7R40, _7R50, _7R60, _7R70; 83 static LUTf _5R10, _5R20, _5R30; 84 static LUTf _25R10, _25R20, _25R30; 85 static LUTf _10RP10, _10RP20, _10RP30; 86 static LUTf _7G30, _7G40, _7G50, _7G60, _7G70, _7G80; 87 static LUTf _5G30, _5G40, _5G50, _5G60, _5G70, _5G80; 88 static LUTf _25G30, _25G40, _25G50, _25G60, _25G70, _25G80; 89 static LUTf _1G30, _1G40, _1G50, _1G60, _1G70, _1G80; 90 static LUTf _10GY30, _10GY40, _10GY50, _10GY60, _10GY70, _10GY80; 91 static LUTf _75GY30, _75GY40, _75GY50, _75GY60, _75GY70, _75GY80; 92 static LUTf _5GY30, _5GY40, _5GY50, _5GY60, _5GY70, _5GY80; 93 94 // Separated from init() to keep the code clear 95 static void initMunsell (); 96 static double hue2rgb(double p, double q, double t); 97 static float hue2rgbfloat(float p, float q, float t); 98 #ifdef __SSE2__ 99 static vfloat hue2rgb(vfloat p, vfloat q, vfloat t); 100 #endif 101 102 static float computeXYZ2Lab(float f); 103 static float computeXYZ2LabY(float f); 104 105 public: 106 107 typedef enum Channel { 108 CHANNEL_RED = 1 << 0, 109 CHANNEL_GREEN = 1 << 1, 110 CHANNEL_BLUE = 1 << 2, 111 CHANNEL_HUE = 1 << 3, 112 CHANNEL_SATURATION = 1 << 4, 113 CHANNEL_VALUE = 1 << 5, 114 CHANNEL_LIGHTNESS = 1 << 6, 115 CHANNEL_CHROMATICITY = 1 << 7 116 } eChannel; 117 118 typedef enum InterpolationPath { 119 IP_SHORTEST, /// Interpolate color using the shortest path between 2 hues 120 IP_LONGEST, /// Interpolate color using the longest path between 2 hues 121 } eInterpolationPath; 122 123 typedef enum InterpolationDirection { 124 ID_UP, /// Interpolate color by increasing the hue value, crossing the upper limit 125 ID_DOWN /// Interpolate color by decreasing the hue value, crossing the lower limit 126 } eInterpolationDirection; 127 128 // Wikipedia sRGB: Unlike most other RGB color spaces, the sRGB gamma cannot be expressed as a single numerical value. 129 // The overall gamma is approximately 2.2, consisting of a linear (gamma 1.0) section near black, and a non-linear section elsewhere involving a 2.4 exponent 130 // and a gamma (slope of log output versus log input) changing from 1.0 through about 2.3. 131 constexpr static double sRGBGamma = 2.2; 132 constexpr static double sRGBGammaCurve = 2.4; 133 134 constexpr static double eps = 216.0 / 24389.0; //0.008856 135 constexpr static double eps_max = MAXVALF * eps; //580.40756; 136 constexpr static double kappa = 24389.0 / 27.0; //903.29630; 137 constexpr static double kappaInv = 27.0 / 24389.0; 138 constexpr static double epsilonExpInv3 = 6.0 / 29.0; 139 140 constexpr static float epsf = eps; 141 constexpr static float kappaf = kappa; 142 constexpr static float kappaInvf = kappaInv; 143 constexpr static float epsilonExpInv3f = epsilonExpInv3; 144 145 constexpr static float D50x = 0.9642f; //0.96422; 146 constexpr static float D50z = 0.8249f; //0.82521; 147 constexpr static double u0 = 4.0 * D50x / (D50x + 15 + 3 * D50z); 148 constexpr static double v0 = 9.0 / (D50x + 15 + 3 * D50z); 149 constexpr static double epskap = 8.0; 150 151 constexpr static float c1By116 = 1.0 / 116.0; 152 constexpr static float c16By116 = 16.0 / 116.0; 153 154 static cmsToneCurve* linearGammaTRC; 155 156 static LUTf cachef; 157 static LUTf cachefy; 158 static LUTf gamma2curve; 159 160 // look-up tables for the standard srgb gamma and its inverse (filled by init()) 161 static LUTf igammatab_srgb; 162 static LUTf igammatab_srgb1; 163 static LUTf gammatab_srgb; 164 static LUTf gammatab_srgb1; 165 166 static LUTf denoiseGammaTab; 167 static LUTf denoiseIGammaTab; 168 169 static LUTf igammatab_24_17; 170 static LUTf gammatab_24_17a; 171 static LUTf gammatab_13_2; 172 static LUTf igammatab_13_2; 173 static LUTf gammatab_115_2; 174 static LUTf igammatab_115_2; 175 static LUTf gammatab_145_3; 176 static LUTf igammatab_145_3; 177 178 // look-up tables for the simple exponential gamma 179 static LUTf gammatab; 180 static LUTuc gammatabThumb; // for thumbnails 181 182 183 static void init (); 184 static void cleanup (); 185 186 187 /** 188 * @brief Extract luminance "sRGB" from red/green/blue values 189 * The range of the r, g and b channel has no importance ([0 ; 1] or [0 ; 65535]...) ; r,g,b can be negatives or > max, but must be in "sRGB" 190 * @param r red channel 191 * @param g green channel 192 * @param b blue channel 193 * @return luminance value 194 */ 195 // xyz_sRGBD65 : conversion matrix from XYZ to sRGB for D65 illuminant: we use diagonal values rgbLuminance(float r,float g,float b)196 static float rgbLuminance(float r, float g, float b) 197 { 198 // WArning: The sum of xyz_sRGBd65[1][] is > 1.0 (i.e. 1.0000001), so we use our own adapted values) 199 // 0.2126729, 0.7151521, 0.0721750 200 return r * 0.2126729f + g * 0.7151521f + b * 0.0721750f; 201 } rgbLuminance(double r,double g,double b)202 static double rgbLuminance(double r, double g, double b) 203 { 204 return r * 0.2126729 + g * 0.7151521 + b * 0.0721750; 205 } 206 rgbLuminance(float r,float g,float b,const double workingspace[3][3])207 static float rgbLuminance(float r, float g, float b, const double workingspace[3][3]) 208 { 209 return r * workingspace[1][0] + g * workingspace[1][1] + b * workingspace[1][2]; 210 } 211 212 #ifdef __SSE2__ rgbLuminance(vfloat r,vfloat g,vfloat b,const vfloat workingspace[3])213 static vfloat rgbLuminance(vfloat r, vfloat g, vfloat b, const vfloat workingspace[3]) 214 { 215 return r * workingspace[0] + g * workingspace[1] + b * workingspace[2]; 216 } 217 #endif 218 219 /** 220 * @brief Convert red/green/blue to L*a*b 221 * @brief Convert red/green/blue to hue/saturation/luminance 222 * @param profile output profile name 223 * @param profileW working profile name 224 * @param r red channel [0 ; 1] 225 * @param g green channel [0 ; 1] 226 * @param b blue channel [0 ; 1] 227 * @param L Lab L channel [0 ; 1] (return value) 228 * @param a Lab a channel [0 ; 1] (return value) 229 * @param b Lab b channel [0; 1] (return value) 230 * @param workingSpace true: compute the Lab value using the Working color space ; false: use the Output color space 231 */ 232 // do not use this function in a loop. It really eats processing time caused by Glib::ustring comparisons 233 static void rgb2lab01 (const Glib::ustring &profile, const Glib::ustring &profileW, float r, float g, float b, float &LAB_l, float &LAB_a, float &LAB_b, bool workingSpace); 234 235 /** 236 * @brief Convert red/green/blue to hue/saturation/luminance 237 * @param r red channel [0 ; 65535] 238 * @param g green channel [0 ; 65535] 239 * @param b blue channel [0 ; 65535] 240 * @param h hue channel [0 ; 1] (return value) 241 * @param s saturation channel [0 ; 1] (return value) 242 * @param l luminance channel [0; 1] (return value) 243 */ 244 static void rgb2hsl (float r, float g, float b, float &h, float &s, float &l); 245 rgb2slfloat(float r,float g,float b,float & s,float & l)246 static inline void rgb2slfloat(float r, float g, float b, float &s, float &l) 247 { 248 249 float minVal = min(r, g, b); 250 float maxVal = max(r, g, b); 251 float C = maxVal - minVal; 252 253 l = (maxVal + minVal) * 7.6295109e-6f; // (0.5f / 65535.f) 254 255 if (C < 0.65535f) { // 0.00001f * 65535.f 256 s = 0.f; 257 } else { 258 259 if (l <= 0.5f) { 260 s = C / (maxVal + minVal); 261 } else { 262 s = C / (131070.f - (maxVal + minVal)); // 131070.f = 2.f * 65535.f 263 } 264 } 265 } 266 rgb2hslfloat(float r,float g,float b,float & h,float & s,float & l)267 static inline void rgb2hslfloat(float r, float g, float b, float &h, float &s, float &l) 268 { 269 270 float minVal = min(r, g, b); 271 float maxVal = max(r, g, b); 272 float C = maxVal - minVal; 273 274 l = (maxVal + minVal) * 7.6295109e-6f; // (0.5f / 65535.f) 275 276 if (C < 0.65535f) { // 0.00001f * 65535.f 277 h = 0.f; 278 s = 0.f; 279 } else { 280 281 if (l <= 0.5f) { 282 s = C / (maxVal + minVal); 283 } else { 284 s = C / (131070.f - (maxVal + minVal)); // 131070.f = 2.f * 65535.f 285 } 286 287 if ( r == maxVal ) { 288 h = (g - b); 289 } else if ( g == maxVal ) { 290 h = (2.f * C) + (b - r); 291 } else { 292 h = (4.f * C) + (r - g); 293 } 294 295 h /= (6.f * C); 296 297 if ( h < 0.f ) { 298 h += 1.f; 299 } 300 } 301 } 302 303 #ifdef __SSE2__ 304 static void rgb2hsl (vfloat r, vfloat g, vfloat b, vfloat &h, vfloat &s, vfloat &l); 305 #endif 306 307 /** 308 * @brief Convert hue/saturation/luminance in red/green/blue 309 * @param h hue channel [0 ; 1] 310 * @param s saturation channel [0 ; 1] 311 * @param l luminance channel [0 ; 1] 312 * @param r red channel [0 ; 65535] (return value) 313 * @param g green channel [0 ; 65535] (return value) 314 * @param b blue channel [0 ; 65535] (return value) 315 */ 316 static void hsl2rgb (float h, float s, float l, float &r, float &g, float &b); 317 hsl2rgbfloat(float h,float s,float l,float & r,float & g,float & b)318 static inline void hsl2rgbfloat (float h, float s, float l, float &r, float &g, float &b) 319 { 320 321 if (s == 0.f) { 322 r = g = b = 65535.f * l; // achromatic 323 } else { 324 float m2; 325 326 if (l <= 0.5f) { 327 m2 = l * (1.f + s); 328 } else { 329 m2 = l + s - l * s; 330 } 331 332 float m1 = 2.f * l - m2; 333 334 r = 65535.f * hue2rgbfloat (m1, m2, h * 6.f + 2.f); 335 g = 65535.f * hue2rgbfloat (m1, m2, h * 6.f); 336 b = 65535.f * hue2rgbfloat (m1, m2, h * 6.f - 2.f); 337 } 338 } 339 340 #ifdef __SSE2__ 341 static void hsl2rgb (vfloat h, vfloat s, vfloat l, vfloat &r, vfloat &g, vfloat &b); 342 #endif 343 344 /** 345 * @brief Convert hue/saturation/luminance in red/green/blue 346 * @param h hue channel [0 ; 1] 347 * @param s saturation channel [0 ; 1] 348 * @param l luminance channel [0 ; 1] 349 * @param r red channel [0 ; 1] (return value) 350 * @param g green channel [0 ; 1] (return value) 351 * @param b blue channel [0 ; 1] (return value) 352 */ 353 static void hsl2rgb01 (float h, float s, float l, float &r, float &g, float &b); 354 355 356 /** 357 * @brief Convert red green blue to hue saturation value 358 * @param r red channel [0 ; 65535] 359 * @param g green channel [0 ; 65535] 360 * @param b blue channel [0 ; 65535] 361 * @param h hue channel [0 ; 1] (return value) 362 * @param s saturation channel [0 ; 1] (return value) 363 * @param v value channel [0 ; 1] (return value) 364 */ 365 static void rgb2hsv (float r, float g, float b, float &h, float &s, float &v); 366 367 /** 368 * @brief Convert red green blue to hue saturation value 369 * @param r red channel [0 ; 1] 370 * @param g green channel [0 ; 1] 371 * @param b blue channel [0 ; 1] 372 * @param h hue channel [0 ; 1] (return value) 373 * @param s saturation channel [0 ; 1] (return value) 374 * @param v value channel [0 ; 1] (return value) 375 */ 376 static void rgb2hsv01 (float r, float g, float b, float &h, float &s, float &v); 377 rgb2s(float r,float g,float b)378 static inline float rgb2s(float r, float g, float b) // fast version if only saturation is needed 379 { 380 float var_Min = min(r, g, b); 381 float var_Max = max(r, g, b); 382 float del_Max = var_Max - var_Min; 383 384 return del_Max / (var_Max == 0.f ? 1.f : var_Max); 385 } 386 rgb2hsvdcp(float r,float g,float b,float & h,float & s,float & v)387 static inline bool rgb2hsvdcp(float r, float g, float b, float &h, float &s, float &v) 388 { 389 390 float var_Min = min(r, g, b); 391 392 if(var_Min < 0.f) { 393 return false; 394 } else { 395 float var_Max = max(r, g, b); 396 float del_Max = var_Max - var_Min; 397 v = var_Max / 65535.f; 398 399 if (fabsf(del_Max) < 0.00001f) { 400 h = 0.f; 401 s = 0.f; 402 } else { 403 s = del_Max / var_Max; 404 405 if ( r == var_Max ) { 406 h = (g - b) / del_Max; 407 } else if ( g == var_Max ) { 408 h = 2.f + (b - r) / del_Max; 409 } else { /*if ( b == var_Max ) */ 410 h = 4.f + (r - g) / del_Max; 411 } 412 413 if ( h < 0.f ) { 414 h += 6.f; 415 } else if ( h > 6.f ) { 416 h -= 6.f; 417 } 418 } 419 420 return true; 421 } 422 } 423 rgb2hsvtc(float r,float g,float b,float & h,float & s,float & v)424 static inline void rgb2hsvtc(float r, float g, float b, float &h, float &s, float &v) 425 { 426 const float var_Min = min(r, g, b); 427 const float var_Max = max(r, g, b); 428 const float del_Max = var_Max - var_Min; 429 430 v = var_Max / 65535.f; 431 432 if (del_Max < 0.00001f) { 433 h = 0.f; 434 s = 0.f; 435 } else { 436 s = del_Max / var_Max; 437 438 if (r == var_Max) { 439 h = (g < b ? 6.f : 0.f) + (g - b) / del_Max; 440 } else if (g == var_Max) { 441 h = 2.f + (b - r) / del_Max; 442 } else { /*if ( b == var_Max ) */ 443 h = 4.f + (r - g) / del_Max; 444 } 445 } 446 } 447 448 /** 449 * @brief Convert hue saturation value in red green blue 450 * @param h hue channel [0 ; 1] 451 * @param s saturation channel [0 ; 1] 452 * @param v value channel [0 ; 1] 453 * @param r red channel [0 ; 65535] (return value) 454 * @param g green channel [0 ; 65535] (return value) 455 * @param b blue channel [0 ; 65535] (return value) 456 */ 457 static void hsv2rgb (float h, float s, float v, float &r, float &g, float &b); 458 hsv2rgbdcp(float h,float s,float v,float & r,float & g,float & b)459 static inline void hsv2rgbdcp (float h, float s, float v, float &r, float &g, float &b) 460 { 461 // special version for dcp which saves 1 division (in caller) and six multiplications (inside this function) 462 const int sector = h; // sector 0 to 5, floor() is very slow, and h is always > 0 463 const float f = h - sector; // fractional part of h 464 465 v *= 65535.f; 466 const float vs = v * s; 467 const float p = v - vs; 468 const float q = v - f * vs; 469 const float t = p + v - q; 470 471 switch (sector) { 472 case 1: 473 r = q; 474 g = v; 475 b = p; 476 break; 477 478 case 2: 479 r = p; 480 g = v; 481 b = t; 482 break; 483 484 case 3: 485 r = p; 486 g = q; 487 b = v; 488 break; 489 490 case 4: 491 r = t; 492 g = p; 493 b = v; 494 break; 495 496 case 5: 497 r = v; 498 g = p; 499 b = q; 500 break; 501 502 default: 503 r = v; 504 g = t; 505 b = p; 506 } 507 } 508 509 static void hsv2rgb (float h, float s, float v, int &r, int &g, int &b); 510 511 512 /** 513 * @brief Convert hue saturation value in red green blue 514 * @param h hue channel [0 ; 1] 515 * @param s saturation channel [0 ; 1] 516 * @param v value channel [0 ; 1] 517 * @param r red channel [0 ; 1] (return value) 518 * @param g green channel [0 ; 1] (return value) 519 * @param b blue channel [0 ; 1] (return value) 520 */ 521 static void hsv2rgb01 (float h, float s, float v, float &r, float &g, float &b); 522 523 524 /** 525 * @brief Convert xyz to red/green/blue 526 * Color space : sRGB - illuminant D50 - use matrix sRGB_xyz[] 527 * @param x X coordinate [0 ; 1] or [0 ; 65535] 528 * @param y Y coordinate [0 ; 1] or [0 ; 65535] 529 * @param z Z coordinate [0 ; 1] or [0 ; 65535] 530 * @param r red channel [same range than xyz channel] (return value) 531 * @param g green channel [same range than xyz channel] (return value) 532 * @param b blue channel [same range than xyz channel] (return value) 533 */ 534 static void xyz2srgb (float x, float y, float z, float &r, float &g, float &b); 535 536 537 /** 538 * @brief Convert xyz to red/green/blue 539 * Color space : Prophoto - illuminant D50 - use the Prophoto_xyz[] matrix 540 * @param x X coordinate [0 ; 1] or [0 ; 65535] 541 * @param y Y coordinate [0 ; 1] or [0 ; 65535] 542 * @param z Z coordinate [0 ; 1] or [0 ; 65535] 543 * @param r red channel [same range than xyz channel] (return value) 544 * @param g green channel [same range than xyz channel] (return value) 545 * @param b blue channel [same range than xyz channel] (return value) 546 */ 547 static void xyz2Prophoto (float x, float y, float z, float &r, float &g, float &b); 548 549 550 /** 551 * @brief Convert rgb in xyz 552 * Color space : Prophoto - illuminant D50 - use matrix xyz_prophoto[] 553 * @param r red channel [0 ; 1] or [0 ; 65535] (return value) 554 * @param g green channel [0 ; 1] or [0 ; 65535] (return value) 555 * @param b blue channel [0 ; 1] or [0 ; 65535] (return value) 556 * @param x X coordinate [same range than xyz channel] 557 * @param y Y coordinate [same range than xyz channel] 558 * @param z Z coordinate [same range than xyz channel] 559 */ 560 static void Prophotoxyz (float r, float g, float b, float &x, float &y, float &z); 561 562 563 /** 564 * @brief Convert xyz in rgb 565 * Color space : undefined - use adhoc matrix: rgb_xyz[3][3] (iccmatrice.h) in function of working space 566 * @param x X coordinate [0 ; 1] or [0 ; 65535] 567 * @param y Y coordinate [0 ; 1] or [0 ; 65535] 568 * @param z Z coordinate [0 ; 1] or [0 ; 65535] 569 * @param r red channel [same range than xyz channel] (return value) 570 * @param g green channel [same range than xyz channel] (return value) 571 * @param b blue channel [same range than xyz channel] (return value) 572 * @param rgb_xyz[3][3] transformation matrix to use for the conversion 573 */ 574 static void xyz2rgb (float x, float y, float z, float &r, float &g, float &b, const double rgb_xyz[3][3]); 575 static void xyz2r (float x, float y, float z, float &r, const double rgb_xyz[3][3]); 576 static void xyz2rgb (float x, float y, float z, float &r, float &g, float &b, const float rgb_xyz[3][3]); 577 #ifdef __SSE2__ 578 static void xyz2rgb (vfloat x, vfloat y, vfloat z, vfloat &r, vfloat &g, vfloat &b, const vfloat rgb_xyz[3][3]); 579 #endif 580 581 582 /** 583 * @brief Convert rgb in xyz 584 * Color space : undefined - use adhoc matrix : xyz_rgb[3][3] (iccmatrice.h) in function of working space 585 * @param r red channel [0 ; 1] or [0 ; 65535] 586 * @param g green channel [0 ; 1] or [0 ; 65535] 587 * @param b blue channel [0 ; 1] or [0 ; 65535] 588 * @param x X coordinate [same range than rgb channel] (return value) 589 * @param y Y coordinate [same range than rgb channel] (return value) 590 * @param z Z coordinate [same range than rgb channel] (return value) 591 * @param xyz_rgb[3][3] transformation matrix to use for the conversion 592 */ 593 static void rgbxyz (float r, float g, float b, float &x, float &y, float &z, const double xyz_rgb[3][3]); 594 static void rgbxyz (float r, float g, float b, float &x, float &y, float &z, const float xyz_rgb[3][3]); 595 #ifdef __SSE2__ 596 static void rgbxyz (vfloat r, vfloat g, vfloat b, vfloat &x, vfloat &y, vfloat &z, const vfloat xyz_rgb[3][3]); 597 #endif 598 599 /** 600 * @brief Convert Lab in xyz 601 * @param L L channel [0 ; 32768] ; L can be negative rarely or superior 32768 602 * @param a channel [-42000 ; +42000] ; can be more than 42000 603 * @param b channel [-42000 ; +42000] ; can be more than 42000 604 * @param x X coordinate [0 ; 65535] ; can be negative! (return value) 605 * @param y Y coordinate [0 ; 65535] ; can be negative! (return value) 606 * @param z Z coordinate [0 ; 65535] ; can be negative! (return value) 607 */ 608 static void Lab2XYZ(float L, float a, float b, float &x, float &y, float &z); 609 static void L2XYZ(float L, float &x, float &y, float &z); 610 611 #ifdef __SSE2__ 612 static void Lab2XYZ(vfloat L, vfloat a, vfloat b, vfloat &x, vfloat &y, vfloat &z); 613 #endif // __SSE2__ 614 615 /** 616 * @brief Convert xyz in Lab 617 * @param x X coordinate [0 ; 65535] ; can be negative or superior to 65535 618 * @param y Y coordinate [0 ; 65535] ; can be negative or superior to 65535 619 * @param z Z coordinate [0 ; 65535] ; can be negative or superior to 65535 620 * @param L L channel [0 ; 32768] ; L can be negative rarely or superior 32768 (return value) 621 * @param a channel [-42000 ; +42000] ; can be more than 42000 (return value) 622 * @param b channel [-42000 ; +42000] ; can be more than 42000 (return value) 623 */ 624 static void XYZ2Lab(float x, float y, float z, float &L, float &a, float &b); 625 static void RGB2Lab(float *X, float *Y, float *Z, float *L, float *a, float *b, const float wp[3][3], int width); 626 static void Lab2RGBLimit(float *L, float *a, float *b, float *R, float *G, float *B, const float wp[3][3], float limit, float afactor, float bfactor, int width); 627 static void RGB2L(float *X, float *Y, float *Z, float *L, const float wp[3][3], int width); 628 629 /** 630 * @brief Convert Lab in Yuv 631 * @param L L channel [0 ; 32768] ; L can be negative rarely or superior 32768 632 * @param a channel [-42000 ; +42000] ; can be more than 42000 633 * @param b channel [-42000 ; +42000] ; can be more than 42000 634 * @param Y luminance channel [0 ; 65535] (return value) 635 * @param u red chrominance channel [0 ; 65535] (return value) 636 * @param v blue chrominance channel [0 ; 65535] (return value) 637 */ 638 static void Lab2Yuv(float L, float a, float b, float &Y, float &u, float &v); 639 640 641 /** 642 * @brief Convert Yuv in Lab 643 * @param Y luminance channel [0 ; 65535] 644 * @param u red chrominance channel [0 ; 65535] 645 * @param v blue chrominance channel [0 ; 65535] 646 * @param L L channel [0 ; 32768] ; L can be negative rarely or superior 32768 (return value) 647 * @param a channel [-42000 ; +42000] ; can be more than 42000 (return value) 648 * @param b channel [-42000 ; +42000] ; can be more than 42000 (return value) 649 */ 650 static void Yuv2Lab(float Y, float u, float v, float &L, float &a, float &b, const double wp[3][3]); 651 652 653 /** 654 * @brief Convert the 'a' and 'b' channels of the L*a*b color space to 'c' and 'h' channels of the Lch color space (channel 'L' is identical [0 ; 32768]) 655 * @param a 'a' channel [-42000 ; +42000] ; can be more than 42000 656 * @param b 'b' channel [-42000 ; +42000] ; can be more than 42000 657 * @param c 'c' channel return value, in [0 ; 42000] ; can be more than 42000 (return value) 658 * @param h 'h' channel return value, in [-PI ; +PI] (return value) 659 */ 660 static void Lab2Lch(float a, float b, float &c, float &h); 661 #ifdef __SSE2__ 662 static void Lab2Lch(float *a, float *b, float *c, float *h, int w); 663 #endif 664 665 /** 666 * @brief Convert 'c' and 'h' channels of the Lch color space to the 'a' and 'b' channels of the L*a*b color space (channel 'L' is identical [0 ; 32768]) 667 * @param c 'c' channel value, in [0 ; 42000] 668 * @param h 'h' channel value, in [-PI ; +PI] 669 * @param a 'a' channel [-42000 ; +42000] ; can be more than 42000 (return value) 670 * @param b 'b' channel [-42000 ; +42000] ; can be more than 42000 (return value) 671 */ 672 static void Lch2Lab(float c, float h, float &a, float &b); 673 674 675 /** 676 * @brief Convert the 'u' and 'v' channels of the Luv color space to 'c' and 'h' channels of the Lch color space ('L' channel is identical) 677 * @param u 'u' channel [unknown range!] 678 * @param v 'v' channel [unknown range!] 679 * @param c 'c' channel [unknown range!] (return value) 680 * @param h 'h' channel [-PI ; +PI] (return value) 681 */ 682 static void Luv2Lch(float u, float v, float &c, float &h); 683 684 685 /** 686 * @brief Convert 'c' and 'h' channels of the Lch color space to the 'u' and 'v' channels of the Luv color space ('L' channel is identical) 687 * @param c 'c' channel [unknown range!] ; can be more than 42000 688 * @param h 'h' channel [-PI ; +PI] 689 * @param u 'u' channel [unknown range!] (return value) 690 * @param v 'v' channel [unknown range!] (return value) 691 */ 692 static void Lch2Luv(float c, float h, float &u, float &v); 693 694 695 /** 696 * @brief Return "f" in function of CIE's kappa and epsilon constants 697 * @param f f can be fx fy fz where: 698 * fx=a/500 + fy a=chroma green red [-128 ; +128] 699 * fy=(L+16)/116 L=luminance [0 ; 100] 700 * fz=fy-b/200 b=chroma blue yellow [-128 ; +128] 701 */ f2xyz(double f)702 static inline double f2xyz(double f) 703 { 704 return (f > epsilonExpInv3) ? f * f * f : (116. * f - 16.) * kappaInv; 705 706 } f2xyz(float f)707 static inline float f2xyz(float f) 708 { 709 return (f > epsilonExpInv3f) ? f * f * f : (116.f * f - 16.f) * kappaInvf; 710 } 711 #ifdef __SSE2__ f2xyz(vfloat f)712 static inline vfloat f2xyz(vfloat f) 713 { 714 const vfloat epsilonExpInv3v = F2V(epsilonExpInv3f); 715 const vfloat kappaInvv = F2V(kappaInvf); 716 vfloat res1 = f * f * f; 717 vfloat res2 = (F2V(116.f) * f - F2V(16.f)) * kappaInvv; 718 return vself(vmaskf_gt(f, epsilonExpInv3v), res1, res2); 719 } 720 #endif 721 722 /** 723 * @brief Calculate the effective direction (up or down) to linearly interpolating 2 colors so that it follows the shortest or longest path 724 * @param h1 First hue [0 ; 1] 725 * @param h2 Second hue [0 ; 1] 726 * @param path Path to follow (shortest/longest) 727 * @return The interpolation direction 728 */ getHueInterpolationDirection(double h1,double h2,eInterpolationPath path)729 static inline eInterpolationDirection getHueInterpolationDirection (double h1, double h2, eInterpolationPath path) 730 { 731 if (path == IP_SHORTEST) { 732 if (h2 > h1) { 733 if (h2 - h1 <= 0.5) { 734 return ID_UP; 735 } else { 736 return ID_DOWN; 737 } 738 } else { 739 if (h1 - h2 <= 0.5) { 740 return ID_DOWN; 741 } else { 742 return ID_UP; 743 } 744 } 745 } else { 746 if (h2 > h1) { 747 if (h2 - h1 <= 0.5) { 748 return ID_DOWN; 749 } else { 750 return ID_UP; 751 } 752 } else { 753 if (h1 - h2 <= 0.5) { 754 return ID_UP; 755 } else { 756 return ID_DOWN; 757 } 758 } 759 } 760 } 761 762 763 /** 764 * @brief Calculate a color by linearly interpolating 2 colors 765 * @param h1 First hue 766 * @param h2 Second hue 767 * @param balance Factor from 0 (first hue) to 1 (second hue) 768 * @param dir Tells which direction the interpolation have to follow. You can get the value with getHueInterpolationDirection 769 * @return The interpolated hue 770 */ interpolateHueHSV(double h1,double h2,double balance,eInterpolationDirection dir)771 static inline double interpolateHueHSV (double h1, double h2, double balance, eInterpolationDirection dir) 772 { 773 if (h1 == h2) { 774 return h1; 775 } 776 777 if (dir == ID_DOWN) { 778 if (h1 < h2) { 779 double temp = h1; 780 h1 = h2 - 1.; 781 h2 = temp; 782 balance = 1. - balance; 783 } 784 785 double h3 = h1 + balance * (h2 - h1); 786 787 if (h3 < 0.) { 788 h3 += 1.; 789 } 790 791 return h3; 792 } else { 793 if (h1 > h2) { 794 h2 += 1.; 795 } 796 797 double h3 = h1 + balance * (h2 - h1); 798 799 if (h3 > 1.) { 800 h3 -= 1.; 801 } 802 803 return h3; 804 } 805 } 806 807 /** 808 * @brief Interpolate 2 colors from their respective red/green/blue channels, with a balance factor 809 * @param balance gives weight to the first and second color [0 ; 1] 810 * 0. = output color == first color 811 * 0.5 = output color == equally mixed colors 812 * 1. = output color == second color 813 * @param r1 red channel of color 1 [0 ; 65535] 814 * @param g1 green channel of color 1 [0 ; 65535] 815 * @param b1 blue channel of color 1 [0 ; 65535] 816 * @param r2 red channel of color 2 [0 ; 65535] 817 * @param g2 green channel of color 2 [0 ; 65535] 818 * @param b2 blue channel of color 2 [0 ; 65535] 819 * @param channels bitfield of channel to interpolate (CHANNEL_LIGHTNESS|CHANNEL_CHROMATICITY|CHANNEL_HUE) 820 * @param xyz_rgb color space 821 * @param ro red channel of output color [0 ; 65535] (return value) 822 * @param go green channel of output color [0 ; 65535] (return value) 823 * @param bo blue channel of output color [0 ; 65535] (return value) 824 */ 825 static void interpolateRGBColor (float balance, float r1, float g1, float b1, float r2, float g2, float b2, int channels, const double xyz_rgb[3][3], const double rgb_xyz[3][3], float &ro, float &go, float &bo); 826 827 /** 828 * @brief Interpolate 2 colors from their respective red/green/blue channels, with a balance factor 829 * @param realL luminance hsl [0; 1] 830 * @param iplow low luminance for rl [0;1] 831 * @param ihigh high luminance for r2 [0;1] 832 * @param algm algorithm [0;2] 833 * @param balance gives weight to the first and second color [0 ; 1] 834 * 0. = output color == first color 835 * 0.5 = output color == equally mixed colors 836 * 1. = output color == second color 837 * @param twoc 2 colors or 512 int 838 * @param r1 red channel of color 1 [0 ; 65535] 839 * @param g1 green channel of color 1 [0 ; 65535] 840 * @param b1 blue channel of color 1 [0 ; 65535] 841 * @param rl red channel of color low [0 ; 65535] 842 * @param gl green channel of color low [0 ; 65535] 843 * @param bl blue channel of color low [0 ; 65535] 844 845 * @param r2 red channel of color 2 or high[0 ; 65535] 846 * @param g2 green channel of color 2 or high[0 ; 65535] 847 * @param b2 blue channel of color 2 [or high 0 ; 65535] 848 * @param channels bitfield of channel to interpolate (CHANNEL_LIGHTNESS|CHANNEL_CHROMATICITY|CHANNEL_HUE) 849 * @param xyz_rgb color space 850 * @param rgb_xyz inverse color space 851 * @param ro red channel of output color [0 ; 65535] (return value) 852 * @param go green channel of output color [0 ; 65535] (return value) 853 * @param bo blue channel of output color [0 ; 65535] (return value) 854 */ 855 static void interpolateRGBColor (float realL, float iplow, float iphigh, int algm, float balance, int twoc, int metchrom, float chromat, float luma, float r1, float g1, float b1, float xl, float yl, float zl, float x2, float y2, float z2, const double xyz_rgb[3][3], const double rgb_xyz[3][3], float &ro, float &go, float &bo); 856 857 858 /** 859 * @brief Interpolate a hue value as the angle of a polar coordinate with hue in the [0;1] range 860 * Chose the shorter path from hue 1 to hue 2. 861 * @param h1 First hue [0; 1] 862 * @param h2 Second hue [0; 1] 863 * @param balance Interpolation factor [0 ; 1] where 0.=h1, 1.=h2 864 * @return the interpolated value [0;1] 865 */ 866 /*template <typename T, typename U> 867 static inline T interpolatePolarHue_01 (T h1, T h2, U balance) { 868 869 if (h1==h2) 870 return h1; 871 if ((h1 > h2) && (h1-h2 > T(0.5))){ 872 h1 -= T(1.); 873 double value = h1 + T(balance) * (h2-h1); 874 if (value < T(0.)) 875 value += T(1.); 876 return value; 877 } 878 else if (h2-h1 > T(0.5)) { 879 h2 -= T(1.); 880 double value = h1 + T(balance) * (h2-h1); 881 if (value < T(0.)) 882 value += T(1.); 883 return value; 884 } 885 else 886 return h1 + T(balance) * (h2-h1); 887 }*/ 888 889 890 /** 891 * @brief Interpolate a hue value as the angle of a polar coordinate with hue in the [-PI ; +PI] range 892 * Chose the shorter path from hue 1 to hue 2. 893 * @param h1 First hue [-PI ; +PI] 894 * @param h2 Second hue [-PI ; +PI] 895 * @param balance Interpolation factor [0 ; 1] where 0.=h1, 1.=h2 896 * @return the interpolated value [-PI ; +PI] 897 */ 898 /*template <typename T, typename U> 899 static inline T interpolatePolarHue_PI (T h1, T h2, U balance) { 900 if (h1==h2) 901 return h1; 902 if ((h1 > h2) && (h1-h2 > T(rtengine::RT_PI))){ 903 h1 -= T(2*rtengine::RT_PI); 904 T value = h1 + T(balance) * (h2-h1); 905 if (value < T(-rtengine::RT_PI)) 906 value += T(2*rtengine::RT_PI); 907 return value; 908 } 909 else if (h2-h1 > T(rtengine::RT_PI)) { 910 h2 -= T(2*rtengine::RT_PI); 911 T value = h1 + T(balance) * (h2-h1); 912 if (value < T(0)) 913 value += T(2*rtengine::RT_PI); 914 return value; 915 } 916 else 917 return h1 + T(balance) * (h2-h1); 918 }*/ 919 920 921 /** 922 * @brief Interpolate a hue value as the angle of a polar coordinate with hue in the [0;1] range 923 * Chose the shorter path from hue 1 to hue 2. 924 * @param h1 First hue [0; 1] 925 * @param h2 Second hue [0; 1] 926 * @param balance Interpolation factor [0 ; 1] where 0.=h1, 1.=h2 927 * @return the interpolated value [0;1] 928 */ 929 template <typename T, typename U> interpolatePolarHue_01(T h1,T h2,U balance)930 static inline T interpolatePolarHue_01 (T h1, T h2, U balance) 931 { 932 float d = h2 - h1; 933 float f; 934 f = T(balance); 935 double h; 936 937 if (h1 > h2) { 938 std::swap(h1, h2); 939 d = -d; 940 f = 1.f - f; 941 } 942 943 if (d < T(-rtengine::RT_PI) || d < T(0) || d > T(rtengine::RT_PI)) { //there was an inversion here !! d > T(rtengine::RT_PI) 944 h1 += T(2 * rtengine::RT_PI); 945 h = h1 + f * (h2 - h1); 946 h = std::fmod(h, 2 * rtengine::RT_PI); 947 } else { 948 h = h1 + f * d; 949 } 950 951 // not strictly necessary..but in case of 952 if(h < T(-rtengine::RT_PI)) { 953 h = T(2 * rtengine::RT_PI) - h; 954 } 955 956 if(h > T(rtengine::RT_PI)) { 957 h = h - T(2 * rtengine::RT_PI); 958 } 959 960 return h; 961 } 962 963 964 /** 965 * @brief Interpolate a hue value as the angle of a polar coordinate with hue in the [-PI ; +PI] range 966 * Chose the shorter path from hue 1 to hue 2. 967 * @param h1 First hue [-PI ; +PI] 968 * @param h2 Second hue [-PI ; +PI] 969 * @param balance Interpolation factor [0 ; 1] where 0.=h1, 1.=h2 970 * @return the interpolated value [-PI ; +PI ] 971 */ 972 template <typename T, typename U> interpolatePolarHue_PI(T h1,T h2,U balance)973 static inline T interpolatePolarHue_PI (T h1, T h2, U balance) 974 { 975 float d = h2 - h1; 976 float f; 977 f = T(balance); 978 double h; 979 980 if (h1 > h2) { 981 std::swap(h1, h2); 982 d = -d; 983 f = 1.f - f; 984 } 985 986 if (d < T(0) || d < T(0.5) || d > T(1.)) { //there was an inversion here !! d > T(rtengine::RT_PI) 987 h1 += T(1.); 988 h = h1 + f * (h2 - h1); 989 h = std::fmod(h, 1.); 990 } else { 991 h = h1 + f * d; 992 } 993 994 // not strictly necessary..but in case of 995 if(h < T(0)) { 996 h = T(1.) - h; 997 } 998 999 if(h > T(1)) { 1000 h = h - T(1.); 1001 } 1002 1003 return h; 1004 } 1005 1006 /** 1007 * @brief Get the gamma curves' parameters used by LCMS2 1008 * @param pwr gamma value [>1] 1009 * @param ts slope [0 ; 20] 1010 * @param mode [always 0] 1011 * @imax imax [always 0] 1012 * @param gamma a pointer to an array of 6 double gamma values: 1013 * gamma0 used in ip2Lab2rgb [0 ; 1], usually near 0.5 (return value) 1014 * gamma1 used in ip2Lab2rgb [0 ; 20], can be superior to 20, but it's quite unusual(return value) 1015 * gamma2 used in ip2Lab2rgb [0 ; 1], usually near 0.03(return value) 1016 * gamma3 used in ip2Lab2rgb [0 ; 1], usually near 0.003(return value) 1017 * gamma4 used in ip2Lab2rgb [0 ; 1], usually near 0.03(return value) 1018 * gamma5 used in ip2Lab2rgb [0 ; 1], usually near 0.5 (return value) 1019 */ 1020 static void calcGamma (double pwr, double ts, int mode, GammaValues &gamma); 1021 1022 1023 /** 1024 * @brief Used by Black and White to correct gamma for each channel red, green and blue channel 1025 * @param r red channel input and output value [0 ; 65535] 1026 * @param g green channel input and output value [0 ; 65535] 1027 * @param b blue channel input and output value [0 ; 65535] 1028 * @param gammabwr gamma value for red channel [>0] 1029 * @param gammabwg gamma value for red channel [>0] 1030 * @param gammabwb gamma value for red channel [>0] 1031 */ 1032 static void trcGammaBW (float &r, float &g, float &b, float gammabwr, float gammabwg, float gammabwb); 1033 #ifdef __SSE2__ 1034 static void trcGammaBWRow (float *r, float *g, float *b, int width, float gammabwr, float gammabwg, float gammabwb); 1035 #endif 1036 1037 1038 /** @brief Compute the B&W constants for the Black and White processing and its GUI 1039 * @param setting main mode 1040 * @param filter string of the filter effect to use 1041 * @param algo choice between linear and special for OYCPM colors 1042 * @param mixerRed red channel value of the channel mixer [-100 ; +200] 1043 * @param mixerGreen green channel value of the channel mixer [-100 ; +200] 1044 * @param mixerBlue blue channel value of the channel mixer [-100 ; +200] 1045 * @param mixerOrange orange channel value of the channel mixer [-100 ; +200] 1046 * @param mixerYellow yellow channel value of the channel mixer [-100 ; +200] 1047 * @param mixerCyan cyan channel value of the channel mixer [-100 ; +200] 1048 * @param mixerPurple purple channel value of the channel mixer [-100 ; +200] 1049 * @param mixerMagenta magenta channel value of the channel mixer [-100 ; +200] 1050 * @param autoc automatic mode of the channel mixer 1051 * @param complement adjust complementary channel 1052 * @param kcorec in absolute mode, value to correct the mixer [1 ; 3], usually near 1 (return value) 1053 * @param rrm red channel of the mixer (return value) 1054 * @param ggm green channel of the mixer (return value) 1055 * @param bbm blue channel of the mixer (return value) 1056 */ 1057 static void computeBWMixerConstants (const Glib::ustring &setting, const Glib::ustring &filter, const Glib::ustring &algo, float &filcor, float &mixerRed, float &mixerGreen, 1058 float &mixerBlue, float mixerOrange, float mixerYellow, float mixerCyan, float mixerPurple, float mixerMagenta, 1059 bool autoc, bool complement, float &kcorec, double &rrm, double &ggm, double &bbm); 1060 1061 1062 // standard srgb gamma and its inverse 1063 1064 /** 1065 * @brief sRGB gamma 1066 * See also calcGamma above with the following values: pwr=2.399 ts=12.92310 mode=0.003041 imax=0.055 1067 * @param x red, green or blue channel's value [0 ; 1] 1068 * @return the gamma modified's value [0 ; 1] 1069 */ gamma2(double x)1070 static inline double gamma2 (double x) // g3 1+g4 1071 { 1072 // return x <= 0.003041 ? x * 12.92310 : 1.055 * exp(log(x) / 2.39990) - 0.055;//calculate with calcgamma 1073 //return x <= 0.0031308 ? x * 12.92310 : 1.055 * exp(log(x) / sRGBGammaCurve) - 0.055;//standard discontinuous 1074 //very small differences between the 2 1075 return x <= 0.003040 ? x * 12.92310 : 1.055 * exp(log(x) / sRGBGammaCurve) - 0.055;//continuous 1076 // return x <= 0.003041 ? x * 12.92310 : 1.055011 * exp(log(x) / sRGBGammaCurve) - 0.055011;//continuous 1077 } 1078 1079 1080 /** 1081 * @brief Inverse sRGB gamma 1082 * See also calcGamma above with the following values: pwr=2.3999 ts=12.92310 mode=0.003041 imax=0.055 1083 * @param x red, green or blue channel's value [0 ; 1] 1084 * @return the inverse gamma modified's value [0 ; 1] 1085 */ igamma2(double x)1086 static inline double igamma2 (double x) //g2 1087 { 1088 // return x <= 0.039289 ? x / 12.92310 : exp(log((x + 0.055) / 1.055) * 2.39990);//calculate with calcgamma 1089 // return x <= 0.04045 ? x / 12.92310 : exp(log((x + 0.055) / 1.055) * sRGBGammaCurve);//standard discontinuous 1090 //very small differences between the 4 1091 return x <= 0.039286 ? x / 12.92310 : exp(log((x + 0.055) / 1.055) * sRGBGammaCurve);//continuous 1092 // return x <= 0.039293 ? x / 12.92310 : exp(log((x + 0.055011) / 1.055011) * sRGBGammaCurve);//continuous 1093 } 1094 1095 1096 /** 1097 * @brief Get the gamma value for Gamma=5.5 Slope=10 1098 * @param x red, green or blue channel's value [0 ; 1] 1099 * @return the gamma modified's value [0 ; 1] 1100 */ gamma55(double x)1101 static inline double gamma55 (double x) // g3 1+g4 1102 { 1103 return x <= 0.013189 ? x * 10.0 : 1.593503 * exp(log(x) / 5.5) - 0.593503; // 5.5 10 1104 } 1105 1106 1107 /** 1108 * @brief Get the inverse gamma value for Gamma=5.5 Slope=10 1109 * @param x red, green or blue channel's value [0 ; 1] 1110 * @return the inverse gamma modified's value [0 ; 1] 1111 */ igamma55(double x)1112 static inline double igamma55 (double x) //g2 1113 { 1114 return x <= 0.131889 ? x / 10.0 : exp(log((x + 0.593503) / 1.593503) * 5.5); // 5.5 10 1115 } 1116 1117 1118 /** 1119 * @brief Get the gamma value for Gamma=4 Slope=5 1120 * @param x red, green or blue channel's value [0 ; 1] 1121 * @return the gamma modified's value [0 ; 1] 1122 */ gamma4(double x)1123 static inline double gamma4 (double x) // g3 1+g4 1124 { 1125 return x <= 0.03089 ? x * 5.0 : 1.478793 * exp(log(x) / 4.1) - 0.478793; // 4 5 1126 } 1127 1128 1129 /** 1130 * @brief Get the inverse gamma value for Gamma=4 Slope=5 1131 * @param x red, green or blue channel's value [0 ; 1] 1132 * @return the inverse gamma modified's value [0 ; 1] 1133 */ igamma4(double x)1134 static inline double igamma4 (double x) //g2 1135 { 1136 return x <= 0.154449 ? x / 5.0 : exp(log((x + 0.478793) / 1.478793) * 4.1); // 4 5 1137 } 1138 1139 1140 /* 1141 * @brief Get the gamma value for Gamma=2.2 Slope=4.5 1142 * @param x red, green or blue channel's value [0 ; 1] 1143 * @return the gamma modified's value [0 ; 1] 1144 * 1145 static inline double gamma709 (double x) { 1146 return x <= 0.0176 ? x*4.5 : 1.0954*exp(log(x)/2.2)-0.0954; 1147 } 1148 1149 * @brief Get the inverse gamma value for Gamma=2.2 Slope=4.5 1150 * @param x red, green or blue channel's value [0 ; 1] 1151 * @return the inverse gamma modified's value [0 ; 1] 1152 * 1153 static inline double igamma709 (double x) { 1154 return x <= 0.0795 ? x/4.5 : exp(log((x+0.0954)/1.0954)*2.2); 1155 } 1156 */ 1157 1158 1159 1160 /** 1161 * @brief Get the gamma value for Gamma=2.4 Slope=17 1162 * @param x red, green or blue channel's value [0 ; 1] 1163 * @return the gamma modified's value [0 ; 1] 1164 */ gamma24_17(double x)1165 static inline double gamma24_17 (double x) 1166 { 1167 return x <= 0.001867 ? x * 17.0 : 1.044445 * exp(log(x) / 2.4) - 0.044445; 1168 } 1169 1170 1171 /** 1172 * @brief Get the inverse gamma value for Gamma=2.4 Slope=17 1173 * @param x red, green or blue channel's value [0 ; 1] 1174 * @return the inverse gamma modified's value [0 ; 1] 1175 */ igamma24_17(double x)1176 static inline double igamma24_17 (double x) 1177 { 1178 return x <= 0.031746 ? x / 17.0 : exp(log((x + 0.044445) / 1.044445) * 2.4); 1179 } 1180 1181 1182 /** 1183 * @brief Get the gamma value for Gamma=2.6 Slope=11 1184 * @param x red, green or blue channel's value [0 ; 1] 1185 * @return the gamma modified's value [0 ; 1] 1186 */ gamma26_11(double x)1187 static inline double gamma26_11 (double x) 1188 { 1189 return x <= 0.004921 ? x * 11.0 : 1.086603 * exp(log(x) / 2.6) - 0.086603; 1190 } 1191 1192 1193 /** 1194 * @brief Get the inverse gamma value for Gamma=2.6 Slope=11 1195 * @param x red, green or blue channel's value [0 ; 1] 1196 * @return the inverse gamma modified's value [0 ; 1] 1197 */ igamma26_11(double x)1198 static inline double igamma26_11 (double x) 1199 { 1200 return x <= 0.054127 ? x / 11.0 : exp(log((x + 0.086603) / 1.086603) * 2.6); 1201 } 1202 /** 1203 * @brief Get the gamma value for Gamma=1.3 Slope=2 1204 * @param x red, green or blue channel's value [0 ; 1] 1205 * @return the gamma modified's value [0 ; 1] 1206 */ gamma13_2(double x)1207 static inline double gamma13_2 (double x) 1208 { 1209 return x <= 0.016613 ? x * 2.0 : 1.009968 * exp(log(x) / 1.3) - 0.009968; 1210 } 1211 igamma13_2(double x)1212 static inline double igamma13_2 (double x) 1213 { 1214 return x <= 0.033226 ? x / 2.0 : exp(log((x + 0.009968) / 1.009968) * 1.3); 1215 } 1216 gamma115_2(double x)1217 static inline double gamma115_2 (double x) 1218 { 1219 return x <= 0.001692 ? x * 2.0 : 1.000508 * exp(log(x) / 1.15) - 0.000508; 1220 } 1221 igamma115_2(double x)1222 static inline double igamma115_2 (double x) 1223 { 1224 return x <= 0.003384 ? x / 2.0 : exp(log((x + 0.000508) / 1.000508) * 1.15); 1225 } 1226 gamma145_3(double x)1227 static inline double gamma145_3 (double x) 1228 { 1229 return x <= 0.009115 ? x * 3.0 : 1.012305 * exp(log(x) / 1.45) - 0.012305; 1230 } 1231 igamma145_3(double x)1232 static inline double igamma145_3 (double x) 1233 { 1234 return x <= 0.027345 ? x / 3.0 : exp(log((x + 0.012305) / 1.012305) * 1.45); 1235 } 1236 1237 //gamma for Retinex gammareti(double x,double gamma,double start,double slope,double mul,double add)1238 static inline double gammareti (double x, double gamma, double start, double slope, double mul, double add) 1239 { 1240 return (x <= start ? x*slope : exp(log(x) / gamma) * mul - add); 1241 } igammareti(double x,double gamma,double start,double slope,double mul,double add)1242 static inline double igammareti (double x, double gamma, double start, double slope, double mul, double add) 1243 { 1244 return (x <= start * slope ? x / slope : exp(log((x + add) / mul) * gamma) ); 1245 } 1246 1247 1248 1249 // gamma function with adjustable parameters 1250 //same as above with values calculate with Calcgamma above 1251 // X range 0..1 gamma(double x,double gamma,double start,double slope,double mul,double add)1252 static inline double gamma (double x, double gamma, double start, double slope, double mul, double add) 1253 { 1254 return (x <= start ? x*slope : exp(log(x) / gamma) * mul - add); 1255 } 1256 gammaf(float x,float gamma,float start,float slope)1257 static inline float gammaf (float x, float gamma, float start, float slope) 1258 { 1259 return x <= start ? x * slope : xexpf(xlogf(x) / gamma); 1260 } 1261 1262 //fills a LUT of size 65536 using gamma with slope... 1263 static void gammaf2lut (LUTf &gammacurve, float gamma, float start, float slope, float divisor, float factor); 1264 igamma(double x,double gamma,double start,double slope,double mul,double add)1265 static inline double igamma (double x, double gamma, double start, double slope, double mul, double add) 1266 { 1267 return (x <= start * slope ? x / slope : exp(log((x + add) / mul) * gamma) ); 1268 } 1269 1270 1271 /** 1272 * @brief Very basic gamma 1273 * @param x red, green or blue channel's value [0 ; 1] 1274 * @param gamma gamma value [1 ; 5] 1275 * @return the gamma modified's value [0 ; 1] 1276 */ gamman(double x,double gamma)1277 static inline double gamman (double x, double gamma) //standard gamma without slope... 1278 { 1279 return exp(log(x) / gamma); 1280 } 1281 1282 /** 1283 * @brief Very basic gamma 1284 * @param x red, green or blue channel's value [0 ; 1] 1285 * @param gamma gamma value [1 ; 5] 1286 * @return the gamma modified's value [0 ; 1] 1287 */ gammanf(float x,float gamma)1288 static inline float gammanf (float x, float gamma) //standard gamma without slope... 1289 { 1290 return xexpf(xlogf(x) / gamma); 1291 } 1292 //fills a LUT of size 65536 using gamma without slope... 1293 static void gammanf2lut (LUTf &gammacurve, float gamma, float divisor, float factor); 1294 1295 /** 1296 * @brief Very simply inverse gamma 1297 * @param x red, green or blue channel's value [0 ; 1] 1298 * @param gamma gamma value [1 ; 5] 1299 * @return the inverse gamma modified's value [0 ; 1] 1300 */ igamman(double x,double gamma)1301 static inline double igamman (double x, double gamma) //standard inverse gamma without slope... 1302 { 1303 return exp(log(x) * gamma); 1304 } 1305 1306 1307 /** 1308 * @brief Get the gamma value out of look-up tables 1309 * Calculated with gamma function above. e.g. : 1310 * for (int i=0; i<65536; i++) 1311 * gammatab_srgb[i] = (65535.0 * gamma2 (i/65535.0)); 1312 * @param x [0 ; 1] 1313 * @return the gamma modified's value [0 ; 65535] 1314 */ gamma_srgb(char x)1315 static inline float gamma_srgb (char x) 1316 { 1317 return gammatab_srgb[x]; 1318 } gamma(char x)1319 static inline float gamma (char x) 1320 { 1321 return gammatab[x]; 1322 } igamma_srgb(char x)1323 static inline float igamma_srgb (char x) 1324 { 1325 return igammatab_srgb[x]; 1326 } gamma_srgb(int x)1327 static inline float gamma_srgb (int x) 1328 { 1329 return gammatab_srgb[x]; 1330 } gamma(int x)1331 static inline float gamma (int x) 1332 { 1333 return gammatab[x]; 1334 } igamma_srgb(int x)1335 static inline float igamma_srgb (int x) 1336 { 1337 return igammatab_srgb[x]; 1338 } gamma_srgb(float x)1339 static inline float gamma_srgb (float x) 1340 { 1341 return gammatab_srgb[x]; 1342 } gamma_srgbclipped(float x)1343 static inline float gamma_srgbclipped (float x) 1344 { 1345 return gamma2curve[x]; 1346 } gamma(float x)1347 static inline float gamma (float x) 1348 { 1349 return gammatab[x]; 1350 } igamma_srgb(float x)1351 static inline float igamma_srgb (float x) 1352 { 1353 return igammatab_srgb[x]; 1354 } 1355 //static inline float gamma_srgb (double x) { return gammatab_srgb[x]; } 1356 //static inline float gamma (double x) { return gammatab[x]; } 1357 //static inline float igamma_srgb (double x) { return igammatab_srgb[x]; } 1358 1359 1360 1361 // -------------------------------- Jacques's Munsell correction 1362 1363 1364 /** 1365 * @brief Corrects the color (hue) depending on chromaticity and luminance changes 1366 * 1367 * To use in a "for" or "do while" statement. 1368 * 1369 * @param lumaMuns true => luminance correction (for delta L > 10) and chroma correction ; false => only chroma 1370 * @param Lprov1 luminance after [0 ; 100] 1371 * @param Loldd luminance before [0 ; 100] 1372 * @param HH hue before [-PI ; +PI] 1373 * @param Chprov1 chroma after [0 ; 180 (can be superior)] 1374 * @param CC chroma before [0 ; 180] 1375 * @param correctionHueChroma hue correction depending on chromaticity (saturation), in radians [0 ; 0.45] (return value) 1376 * @param correctlum hue correction depending on luminance (brightness, contrast,...), in radians [0 ; 0.45] (return value) 1377 * @param munsDbgInfo (Debug target only) object to collect information 1378 */ 1379 1380 #ifdef _DEBUG 1381 static void AllMunsellLch (bool lumaMuns, float Lprov1, float Loldd, float HH, float Chprov1, float CC, float &correctionHueChroma, float &correctlum, MunsellDebugInfo* munsDbgInfo); 1382 #else 1383 static void AllMunsellLch (bool lumaMuns, float Lprov1, float Loldd, float HH, float Chprov1, float CC, float &correctionHueChroma, float &correctlum); 1384 #endif 1385 static void AllMunsellLch (float Lprov1, float HH, float Chprov1, float CC, float &correctionHueChroma); 1386 1387 1388 /** 1389 * @brief Correct chromaticity and luminance so that the color stays in the working profile's gamut 1390 * 1391 * This function puts the data (Lab) in the gamut of "working profile": 1392 * it returns the corrected values of the chromaticity and luminance 1393 * 1394 * @param HH : hue, in radians [-PI ; +PI] 1395 * @param Lprov1 : input luminance value, sent back corrected [0 ; 100] (input & output value) 1396 * @param Chprov1: input chroma value, sent back corrected [0 ; 180 (can be superior)] (input & output value) 1397 * @param R red value of the corrected color [0 ; 65535 but can be negative or superior to 65535] (return value) 1398 * @param G green value of the corrected color [0 ; 65535 but can be negative or superior to 65535] (return value) 1399 * @param B blue value of the corrected color [0 ; 65535 but can be negative or superior to 65535] (return value) 1400 * @param wip working profile 1401 * @param isHLEnabled true if "Highlight Reconstruction " is enabled 1402 * @param lowerCoef a float number between [0.95 ; 1.0[ 1403 * The nearest it is from 1.0, the more precise it will be, and the longer too as more iteration will be necessary 1404 * @param higherCoef a float number between [0.95 ; 1.0[ 1405 * The nearest it is from 1.0, the more precise it will be, and the longer too as more iteration will be necessary 1406 * @param neg (Debug target only) to calculate iterations for negatives values 1407 * @param moreRGB (Debug target only) to calculate iterations for values >65535 1408 */ 1409 #ifdef _DEBUG 1410 static void gamutLchonly (float HH, float &Lprov1, float &Chprov1, float &R, float &G, float &B, const double wip[3][3], bool isHLEnabled, float lowerCoef, float higherCoef, bool &neg, bool &more_rgb); 1411 static void gamutLchonly (float HH, float2 sincosval, float &Lprov1, float &Chprov1, float &R, float &G, float &B, const double wip[3][3], bool isHLEnabled, float lowerCoef, float higherCoef, bool &neg, bool &more_rgb); 1412 static void gamutLchonly (float2 sincosval, float &Lprov1, float &Chprov1, const float wip[3][3], bool isHLEnabled, float lowerCoef, float higherCoef, bool &neg, bool &more_rgb); 1413 #else 1414 static void gamutLchonly (float HH, float &Lprov1, float &Chprov1, float &R, float &G, float &B, const double wip[3][3], bool isHLEnabled, float lowerCoef, float higherCoef); 1415 static void gamutLchonly (float HH, float2 sincosval, float &Lprov1, float &Chprov1, float &R, float &G, float &B, const double wip[3][3], bool isHLEnabled, float lowerCoef, float higherCoef); 1416 static void gamutLchonly (float2 sincosval, float &Lprov1, float &Chprov1, const float wip[3][3], bool isHLEnabled, float lowerCoef, float higherCoef); 1417 #endif 1418 static void gamutLchonly (float HH, float2 sincosval, float &Lprov1, float &Chprov1, float &saturation, const float wip[3][3], bool isHLEnabled, float lowerCoef, float higherCoef); 1419 1420 1421 /** 1422 * @brief Munsell gamut correction 1423 * 1424 * This function is the overall Munsell's corrections, but only on global statement. It may be better to use local statement with AllMunsellLch. 1425 * They are named accordingly : gamutLchonly and AllMunsellLch 1426 * It can be used before and after treatment (saturation, gamma, luminance, ...) 1427 * 1428 * @param labL L channel input and output image 1429 * L channel's usual range is [0 ; 100], but values can be negative or >100 1430 * @param laba a channel input and output image 1431 * @param labb b channel input and output image 1432 * a and b channel's range is usually [-128 ; +128], but values can be >128 1433 * @param N Number of pixels to process 1434 * @param corMunsell performs Munsell correction 1435 * @param lumaMuns whether to apply luma correction or not (used only if corMuns=true) 1436 * true: apply luma + chroma Munsell correction if delta L > 10; 1437 * false: leaves luma untouched 1438 * @param gamut performs gamutLch 1439 * @param wip matrix for working profile 1440 * @param multiThread whether to parallelize the loop or not 1441 */ 1442 static void LabGamutMunsell (float *labL, float *laba, float *labb, int N, bool corMunsell, bool lumaMuns, bool isHLEnabled, bool gamut, const double wip[3][3]); 1443 1444 1445 /* 1446 * @brief Skin tone protection factor 1447 * Skin colors: mixed from NX2 skin color palette, Von Luschan, and photos of white, black, yellow people... 1448 * There are some little exceptions, but it should cover 99% case. 1449 * Pay attention to white balance, and do not change hue and saturation, upstream of the modification 1450 * Used by vibrance 1451 * @param lum luma value [0 ; 100] 1452 * @param hue hue value [-PI ; +PI] 1453 * @param chrom chroma value [0 ; 180] 1454 * @param satreduc [0.1 ; 1] (return value) 1455 * @param chromx [0 or 1], actually only 0 is used 1456 */ 1457 static void SkinSat (float lum, float hue, float chrom, float &satreduc);//jacques Skin color 1458 1459 1460 /** 1461 * @brief Munsell Lch correction 1462 * Find the right LUT and calculate the correction 1463 * @param lum luma value [0 ; 100] 1464 * @param hue hue value [-PI ; +PI] 1465 * @param chrom chroma value [0 ; 180] 1466 * @param memChprov store chroma [0 ; 180] 1467 * @param correction correction value, in radians [0 ; 0.45] 1468 * @param lbe hue in function of chroma, in radian [-PI ; +PI] 1469 * @param zone [1 ; 4] 1=PB correction + sky 2=red yellow correction 3=Green yellow correction 4=Red purple correction 1470 * @param correctL true=enable the Luminance correction 1471 */ 1472 static void MunsellLch (float lum, float hue, float chrom, float memChprov, float &correction, int zone, float &lbe, bool &correctL);//jacques: Munsell correction 1473 1474 1475 // -------------------------------- end Munsell 1476 1477 1478 static void scalered ( float rstprotection, float param, float limit, float HH, float deltaHH, float &scale, float &scaleext); 1479 static void transitred (float HH, float Chprov1, float dred, float factorskin, float protect_red, float factorskinext, float deltaHH, float factorsat, float &factor); 1480 static void skinred ( double J, double h, double sres, double Sp, float dred, float protect_red, int sk, float rstprotection, float ko, double &s); 1481 static void skinredfloat ( float J, float h, float sres, float Sp, float dred, float protect_red, int sk, float rstprotection, float ko, float &s); 1482 // static void scaleredcdbl ( float skinprot, float param, float limit, float HH, float deltaHH, float &scale,float &scaleext); 1483 SkinSatCbdl(float lum,float hue,float chrom,float skinprot,float & scale,bool neg,float b_l,float t_l,float t_r)1484 static inline void SkinSatCbdl (float lum, float hue, float chrom, float skinprot, float &scale, bool neg, float b_l, float t_l, float t_r) 1485 { 1486 1487 static const float C9 = 8.f, C8 = 15.f, C7 = 12.f, C4 = 7.f, C3 = 5.f, C2 = 5.f, C1 = 5.f; 1488 static const float H9 = 0.05f, H8 = 0.25f, H7 = 0.1f, H4 = 0.02f, H3 = 0.02f, H2 = 0.1f, H1 = 0.1f, H10 = -0.2f, H11 = -0.2f; 1489 1490 // "real" skin color : take into account a slight usage of contrast and saturation in RT if option "skin" = 1, uses implicit factor 1.0 1491 // wide area skin color, useful if not accurate colorimetry or if the user has changed hue and saturation, uses explicit factor 0.6 1492 // wide area for transition, uses explicit factor 0.4 1493 1494 if (lum >= 85.0f) { 1495 if((hue > (t_l + 0.53f - H9) && hue < (t_r + H9)) && (chrom > 8.0f && chrom < (14.0f + C9))) { 1496 scale = (100.f - skinprot) / 100.1f; 1497 } else if (lum >= 92.0f) { 1498 if((hue > t_l + 0.4f && hue < t_r) && (chrom > 7.0f && chrom < (15.0f))) { 1499 scale = (100.f - skinprot * 0.6f) / 100.1f; 1500 } else if ((hue > b_l && hue < t_r) && (chrom > 7.0f && chrom < (18.0f))) { 1501 scale = (100.f - skinprot * 0.4f) / 100.1f; 1502 } 1503 } else if ((hue > t_l + 0.4f && hue < t_r - 0.3f) && (chrom > 7.0f && chrom < (26.0f + C9))) { 1504 scale = (100.f - skinprot * 0.6f) / 100.1f; 1505 } else if ((hue > b_l + 0.05f && hue < t_r) && (chrom > 7.0f && chrom < (35.0f + C9))) { 1506 scale = (100.f - skinprot * 0.4f) / 100.1f; 1507 } 1508 } else if (lum >= 70.0f) { 1509 if((hue > t_l + 0.15f && hue < (t_r - 0.2f + H8)) && (chrom > 8.0f && chrom < (35.0f + C8))) { 1510 scale = (100.f - skinprot) / 100.1f; 1511 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r - 0.2f) && (chrom > 7.0f && chrom < (48.0f + C9) )) { 1512 scale = (100.f - skinprot * 0.6f) / 100.1f; 1513 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r) && (chrom > 7.0f && chrom < (55.0f + C9) )) { 1514 scale = (100.f - skinprot * 0.4f) / 100.1f; 1515 } 1516 } else if (lum >= 52.0f) { 1517 if((hue > t_l && hue < (t_r + H7)) && (chrom > 11.0f && chrom < (35.0f + C7))) { 1518 scale = (100.f - skinprot) / 100.1f; 1519 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r - 0.2f) && (chrom > 7.0f && chrom < (48.0f + C9) )) { 1520 scale = (100.f - skinprot * 0.6f) / 100.1f; 1521 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r) && (chrom > 7.0f && chrom < (55.0f + C9) )) { 1522 scale = (100.f - skinprot * 0.4f) / 100.1f; 1523 } 1524 } else if (lum >= 35.0f) { 1525 if((hue > t_l && hue < (t_r + H4)) && (chrom > 13.0f && chrom < (37.0f + C4))) { 1526 scale = (100.f - skinprot) / 100.1f; 1527 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r - 0.2f) && (chrom > 7.0f && chrom < (48.0f + C9) )) { 1528 scale = (100.f - skinprot * 0.6f) / 100.1f; 1529 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r) && (chrom > 7.0f && chrom < (55.0f + C9) )) { 1530 scale = (100.f - skinprot * 0.4f) / 100.1f; 1531 } 1532 } else if (lum >= 20.0f) { 1533 if((hue > t_l && hue < (t_r + H3)) && (chrom > 7.0f && chrom < (35.0f + C3) )) { 1534 scale = (100.f - skinprot) / 100.1f; 1535 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r - 0.2f) && (chrom > 7.0f && chrom < (48.0f + C9) )) { 1536 scale = (100.f - skinprot * 0.6f) / 100.1f; 1537 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r) && (chrom > 7.0f && chrom < (55.0f + C9) )) { 1538 scale = (100.f - skinprot * 0.4f) / 100.1f; 1539 } 1540 } else if (lum >= 10.0f) { 1541 if((hue > (t_l - 0.25f + H10) && hue < (t_r - 0.3f + H2)) && (chrom > 8.0f && chrom < (23.0f + C2))) { 1542 scale = (100.f - skinprot) / 100.1f; 1543 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r - 0.2f) && (chrom > 7.0f && chrom < (35.0f + C1) )) { 1544 scale = (100.f - skinprot * 0.6f) / 100.1f; 1545 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r - 0.1f) && (chrom > 7.0f && chrom < (45.0f + C1) )) { 1546 scale = (100.f - skinprot * 0.4f) / 100.1f; 1547 } 1548 } else if ((hue > (t_l - 0.2f + H10) && hue < (t_r - 0.3f + H1)) && (chrom > 8.0f && chrom < (23.0f + C1))) { 1549 scale = (100.f - skinprot) / 100.1f; 1550 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r - 0.2f) && (chrom > 7.0f && chrom < (35.0f + C1) )) { 1551 scale = (100.f - skinprot * 0.6f) / 100.1f; 1552 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r - 0.1f) && (chrom > 7.0f && chrom < (45.0f + C1) )) { 1553 scale = (100.f - skinprot * 0.4f) / 100.1f; 1554 } 1555 1556 //extended zone for hair, beard and if user adjust high value for skinprot 1557 if(skinprot > 85.f && chrom < 20.f && neg) { 1558 float modula = -0.0666f * skinprot + 6.66f; 1559 scale *= modula; 1560 } 1561 } 1562 SkinSatCbdl2(float lum,float hue,float chrom,float skinprot,float & scale,bool neg,float b_l,float t_l,float t_r,float b_r,int basc)1563 static inline void SkinSatCbdl2 (float lum, float hue, float chrom, float skinprot, float &scale, bool neg, float b_l, float t_l, float t_r, float b_r, int basc) 1564 { 1565 1566 static const float C9 = 8.f, C8 = 15.f, C7 = 12.f, C4 = 7.f, C3 = 5.f, C2 = 5.f, C1 = 5.f; 1567 static const float H9 = 0.05f, H8 = 0.25f, H7 = 0.1f, H4 = 0.02f, H3 = 0.02f, H2 = 0.1f, H1 = 0.1f, H10 = -0.2f, H11 = -0.2f; 1568 1569 // "real" skin color : take into account a slight usage of contrast and saturation in RT if option "skin" = 1, uses implicit factor 1.0 1570 // wide area skin color, useful if not accurate colorimetry or if the user has changed hue and saturation, uses explicit factor 0.6 1571 // wide area for transition, uses explicit factor 0.4 1572 if((b_l > -0.3f && b_r < 2.f) || basc == 0) { //range maxi skin 1573 if (lum >= 85.0f) { 1574 if((hue > (t_l + 0.53f - H9) && hue < (t_r + H9)) && (chrom > 8.0f && chrom < (14.0f + C9))) { 1575 scale = (100.f - skinprot) / 100.1f; 1576 } else if (lum >= 92.0f) { 1577 if((hue > t_l + 0.4f && hue < t_r) && (chrom > 7.0f && chrom < (15.0f))) { 1578 scale = (100.f - skinprot * 0.6f) / 100.1f; 1579 } else if ((hue > b_l && hue < t_r) && (chrom > 7.0f && chrom < (18.0f))) { 1580 scale = (100.f - skinprot * 0.4f) / 100.1f; 1581 } 1582 } else if ((hue > t_l + 0.4f && hue < t_r - 0.3f) && (chrom > 7.0f && chrom < (26.0f + C9))) { 1583 scale = (100.f - skinprot * 0.6f) / 100.1f; 1584 } else if ((hue > b_l + 0.05f && hue < t_r) && (chrom > 7.0f && chrom < (35.0f + C9))) { 1585 scale = (100.f - skinprot * 0.4f) / 100.1f; 1586 } 1587 } else if (lum >= 70.0f) { 1588 if((hue > t_l + 0.15f && hue < (t_r - 0.2f + H8)) && (chrom > 8.0f && chrom < (35.0f + C8))) { 1589 scale = (100.f - skinprot) / 100.1f; 1590 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r - 0.2f) && (chrom > 7.0f && chrom < (48.0f + C9) )) { 1591 scale = (100.f - skinprot * 0.6f) / 100.1f; 1592 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r) && (chrom > 7.0f && chrom < (55.0f + C9) )) { 1593 scale = (100.f - skinprot * 0.4f) / 100.1f; 1594 } 1595 } else if (lum >= 52.0f) { 1596 if((hue > t_l && hue < (t_r + H7)) && (chrom > 11.0f && chrom < (35.0f + C7))) { 1597 scale = (100.f - skinprot) / 100.1f; 1598 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r - 0.2f) && (chrom > 7.0f && chrom < (48.0f + C9) )) { 1599 scale = (100.f - skinprot * 0.6f) / 100.1f; 1600 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r) && (chrom > 7.0f && chrom < (55.0f + C9) )) { 1601 scale = (100.f - skinprot * 0.4f) / 100.1f; 1602 } 1603 } else if (lum >= 35.0f) { 1604 if((hue > t_l && hue < (t_r + H4)) && (chrom > 13.0f && chrom < (37.0f + C4))) { 1605 scale = (100.f - skinprot) / 100.1f; 1606 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r - 0.2f) && (chrom > 7.0f && chrom < (48.0f + C9) )) { 1607 scale = (100.f - skinprot * 0.6f) / 100.1f; 1608 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r) && (chrom > 7.0f && chrom < (55.0f + C9) )) { 1609 scale = (100.f - skinprot * 0.4f) / 100.1f; 1610 } 1611 } else if (lum >= 20.0f) { 1612 if((hue > t_l && hue < (t_r + H3)) && (chrom > 7.0f && chrom < (35.0f + C3) )) { 1613 scale = (100.f - skinprot) / 100.1f; 1614 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r - 0.2f) && (chrom > 7.0f && chrom < (48.0f + C9) )) { 1615 scale = (100.f - skinprot * 0.6f) / 100.1f; 1616 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r) && (chrom > 7.0f && chrom < (55.0f + C9) )) { 1617 scale = (100.f - skinprot * 0.4f) / 100.1f; 1618 } 1619 } else if (lum >= 10.0f) { 1620 if((hue > (t_l - 0.25f + H10) && hue < (t_r - 0.3f + H2)) && (chrom > 8.0f && chrom < (23.0f + C2))) { 1621 scale = (100.f - skinprot) / 100.1f; 1622 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r - 0.2f) && (chrom > 7.0f && chrom < (35.0f + C1) )) { 1623 scale = (100.f - skinprot * 0.6f) / 100.1f; 1624 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r - 0.1f) && (chrom > 7.0f && chrom < (45.0f + C1) )) { 1625 scale = (100.f - skinprot * 0.4f) / 100.1f; 1626 } 1627 } else if ((hue > (t_l - 0.2f + H10) && hue < (t_r - 0.3f + H1)) && (chrom > 8.0f && chrom < (23.0f + C1))) { 1628 scale = (100.f - skinprot) / 100.1f; 1629 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r - 0.2f) && (chrom > 7.0f && chrom < (35.0f + C1) )) { 1630 scale = (100.f - skinprot * 0.6f) / 100.1f; 1631 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r - 0.1f) && (chrom > 7.0f && chrom < (45.0f + C1) )) { 1632 scale = (100.f - skinprot * 0.4f) / 100.1f; 1633 } 1634 1635 //extended zone for hair, beard and if user adjust high value for skinprot 1636 if(skinprot > 85.f && chrom < 20.f && neg) { 1637 float modula = -0.0666f * skinprot + 6.66f; 1638 scale *= modula; 1639 } 1640 } 1641 //end hue skin algo 1642 else if (basc == 1) { //not hue skin linear transition or mod chroma curve 1643 if(hue >= t_l && hue <= t_r) { 1644 scale = (100.f - skinprot) / 100.1f; 1645 } else if(hue > b_l && hue < t_l) { 1646 float sc = (100.f - skinprot) / 100.1f; 1647 float aa = (1.f - sc) / (b_l - t_l); 1648 float bb = 1.f - aa * b_l; 1649 scale = aa * hue + bb; 1650 } else if(hue > t_r && hue < b_r) { 1651 float sc = (100.f - skinprot) / 100.1f; 1652 float aa = (sc - 1.f) / (t_r - b_r); 1653 float bb = 1.f - aa * b_r; 1654 scale = aa * hue + bb; 1655 } 1656 } 1657 } 1658 1659 1660 SkinSatCbdlCam(float lum,float hue,float chrom,float skinprot,float & scale,bool neg,float b_l,float t_l,float t_r)1661 static inline void SkinSatCbdlCam (float lum, float hue, float chrom, float skinprot, float &scale, bool neg, float b_l, float t_l, float t_r) 1662 { 1663 1664 static const float C9 = 8.f, C8 = 15.f, C7 = 12.f, C4 = 7.f, C3 = 5.f, C2 = 5.f, C1 = 5.f; 1665 static const float H9 = 0.05f, H8 = 0.25f, H7 = 0.1f, H4 = 0.02f, H3 = 0.02f, H2 = 0.1f, H1 = 0.1f, H10 = -0.2f, H11 = -0.2f; 1666 1667 float HH = 0.f; 1668 1669 if (hue > 8.6f && hue <= 74.f ) { 1670 HH = (1.15f / 65.4f) * hue - 0.0012f; //H > 0.15 H<1.3 1671 } else if(hue > 0.f && hue <= 8.6f ) { 1672 HH = (0.19f / 8.6f ) * hue - 0.04f; //H>-0.04 H < 0.15 1673 } else if(hue > 355.f && hue <= 360.f) { 1674 HH = (0.11f / 5.0f ) * hue - 7.96f; //H>-0.15 <-0.04 1675 } else if(hue > 74.f && hue < 95.f ) { 1676 HH = (0.30f / 21.0f) * hue + 0.24285f; //H>1.3 H<1.6 1677 } else if(hue >= 95.f && hue < 137.5f) { 1678 HH = 0.01882f * hue - 0.18823f; // H>1.6 H<2.4 1679 } else if(hue > 285.f && hue <= 355.f) { 1680 HH = 0.1642f * hue - 5.982f; //HH>-1.3 HH <-0.15 1681 } 1682 1683 hue = HH; 1684 1685 // "real" skin color : take into account a slight usage of contrast and saturation in RT if option "skin" = 1, uses implicit factor 1.0 1686 // wide area skin color, useful if not accurate colorimetry or if the user has changed hue and saturation, uses explicit factor 0.6 1687 // wide area for transition, uses explicit factor 0.4 1688 1689 if (lum >= 85.0f) { 1690 if((hue > (t_l + 0.53f - H9) && hue < (t_r + H9)) && (chrom > 8.0f && chrom < (14.0f + C9))) { 1691 scale = (100.f - skinprot) / 100.1f; 1692 } else if (lum >= 92.0f) { 1693 if((hue > t_l + 0.4f && hue < t_r) && (chrom > 7.0f && chrom < (15.0f))) { 1694 scale = (100.f - skinprot * 0.6f) / 100.1f; 1695 } else if ((hue > b_l && hue < t_r) && (chrom > 7.0f && chrom < (18.0f))) { 1696 scale = (100.f - skinprot * 0.4f) / 100.1f; 1697 } 1698 } else if ((hue > t_l + 0.4f && hue < t_r - 0.3f) && (chrom > 7.0f && chrom < (26.0f + C9))) { 1699 scale = (100.f - skinprot * 0.6f) / 100.1f; 1700 } else if ((hue > b_l + 0.05f && hue < t_r) && (chrom > 7.0f && chrom < (35.0f + C9))) { 1701 scale = (100.f - skinprot * 0.4f) / 100.1f; 1702 } 1703 } else if (lum >= 70.0f) { 1704 if((hue > t_l + 0.15f && hue < (t_r - 0.2f + H8)) && (chrom > 8.0f && chrom < (35.0f + C8))) { 1705 scale = (100.f - skinprot) / 100.1f; 1706 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r - 0.2f) && (chrom > 7.0f && chrom < (48.0f + C9) )) { 1707 scale = (100.f - skinprot * 0.6f) / 100.1f; 1708 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r) && (chrom > 7.0f && chrom < (55.0f + C9) )) { 1709 scale = (100.f - skinprot * 0.4f) / 100.1f; 1710 } 1711 } else if (lum >= 52.0f) { 1712 if((hue > t_l && hue < (t_r + H7)) && (chrom > 11.0f && chrom < (35.0f + C7))) { 1713 scale = (100.f - skinprot) / 100.1f; 1714 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r - 0.2f) && (chrom > 7.0f && chrom < (48.0f + C9) )) { 1715 scale = (100.f - skinprot * 0.6f) / 100.1f; 1716 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r) && (chrom > 7.0f && chrom < (55.0f + C9) )) { 1717 scale = (100.f - skinprot * 0.4f) / 100.1f; 1718 } 1719 } else if (lum >= 35.0f) { 1720 if((hue > t_l && hue < (t_r + H4)) && (chrom > 13.0f && chrom < (37.0f + C4))) { 1721 scale = (100.f - skinprot) / 100.1f; 1722 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r - 0.2f) && (chrom > 7.0f && chrom < (48.0f + C9) )) { 1723 scale = (100.f - skinprot * 0.6f) / 100.1f; 1724 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r) && (chrom > 7.0f && chrom < (55.0f + C9) )) { 1725 scale = (100.f - skinprot * 0.4f) / 100.1f; 1726 } 1727 } else if (lum >= 20.0f) { 1728 if((hue > t_l && hue < (t_r + H3)) && (chrom > 7.0f && chrom < (35.0f + C3) )) { 1729 scale = (100.f - skinprot) / 100.1f; 1730 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r - 0.2f) && (chrom > 7.0f && chrom < (48.0f + C9) )) { 1731 scale = (100.f - skinprot * 0.6f) / 100.1f; 1732 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r) && (chrom > 7.0f && chrom < (55.0f + C9) )) { 1733 scale = (100.f - skinprot * 0.4f) / 100.1f; 1734 } 1735 } else if (lum >= 10.0f) { 1736 if((hue > (t_l - 0.25f + H10) && hue < (t_r - 0.3f + H2)) && (chrom > 8.0f && chrom < (23.0f + C2))) { 1737 scale = (100.f - skinprot) / 100.1f; 1738 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r - 0.2f) && (chrom > 7.0f && chrom < (35.0f + C1) )) { 1739 scale = (100.f - skinprot * 0.6f) / 100.1f; 1740 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r - 0.1f) && (chrom > 7.0f && chrom < (45.0f + C1) )) { 1741 scale = (100.f - skinprot * 0.4f) / 100.1f; 1742 } 1743 } else if ((hue > (t_l - 0.2f + H10) && hue < (t_r - 0.3f + H1)) && (chrom > 8.0f && chrom < (23.0f + C1))) { 1744 scale = (100.f - skinprot) / 100.1f; 1745 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r - 0.2f) && (chrom > 7.0f && chrom < (35.0f + C1) )) { 1746 scale = (100.f - skinprot * 0.6f) / 100.1f; 1747 } else if ((hue > (b_l + 0.07f + H11) && hue < t_r - 0.1f) && (chrom > 7.0f && chrom < (45.0f + C1) )) { 1748 scale = (100.f - skinprot * 0.4f) / 100.1f; 1749 } 1750 1751 //extended zone for hair, beard and if user adjust high value for skinprot 1752 if(skinprot > 85.f && chrom < 20.f && neg) { 1753 float modula = -0.0666f * skinprot + 6.66f; 1754 scale *= modula; 1755 } 1756 } 1757 1758 1759 /** 1760 * @brief Gamut correction in the XYZ color space 1761 * @param X X channel input value and corrected output value [0 ; 65535] 1762 * @param Y Y channel input value and corrected output value [0 ; 65535] 1763 * @param Z Z channel input value and corrected output value [0 ; 65535] 1764 * @param p working profile 1765 */ 1766 static void gamutmap(float &X, float &Y, float &Z, const double p[3][3]); 1767 1768 1769 /** 1770 * @brief Get HSV's hue from the Lab's hue 1771 * @param HH Lab's hue value, in radians [-PI ; +PI] 1772 * @return HSV's hue value [0 ; 1] 1773 */ huelab_to_huehsv2(float HH)1774 static inline double huelab_to_huehsv2 (float HH) 1775 { 1776 //hr=translate Hue Lab value (-Pi +Pi) in approximative hr (hsv values) (0 1) [red 1/6 yellow 1/6 green 1/6 cyan 1/6 blue 1/6 magenta 1/6 ] 1777 // with multi linear correspondences (I expect there is no error !!) 1778 double hr = 0.0; 1779 //always put h between 0 and 1 1780 1781 if (HH >= 0.f && HH < 0.6f ) { 1782 hr = 0.11666 * double(HH) + 0.93; //hr 0.93 1.00 full red 1783 } else if (HH >= 0.6f && HH < 1.4f ) { 1784 hr = 0.1125 * double(HH) - 0.0675; //hr 0.00 0.09 red yellow orange 1785 } else if (HH >= 1.4f && HH < 2.f ) { 1786 hr = 0.2666 * double(HH) - 0.2833; //hr 0.09 0.25 orange yellow 1787 } else if (HH >= 2.f && HH < 3.14159f) { 1788 hr = 0.1489 * double(HH) - 0.04785; //hr 0.25 0.42 yellow green green 1789 } else if (HH >= -3.14159f && HH < -2.8f ) { 1790 hr = 0.23419 * double(HH) + 1.1557; //hr 0.42 0.50 green 1791 } else if (HH >= -2.8f && HH < -2.3f ) { 1792 hr = 0.16 * double(HH) + 0.948; //hr 0.50 0.58 cyan 1793 } else if (HH >= -2.3f && HH < -0.9f ) { 1794 hr = 0.12143 * double(HH) + 0.85928; //hr 0.58 0.75 blue blue-sky 1795 } else if (HH >= -0.9f && HH < -0.1f ) { 1796 hr = 0.2125 * double(HH) + 0.94125; //hr 0.75 0.92 purple magenta 1797 } else if (HH >= -0.1f && HH < 0.f ) { 1798 hr = 0.1 * double(HH) + 0.93; //hr 0.92 0.93 red 1799 } 1800 1801 // in case of ! 1802 if (hr < 0.0) { 1803 hr += 1.0; 1804 } else if(hr > 1.0) { 1805 hr -= 1.0; 1806 } 1807 1808 return (hr); 1809 } 1810 RGB2Y(const float * R,const float * G,const float * B,float * Y1,float * Y2,int W)1811 static inline void RGB2Y(const float* R, const float* G, const float* B, float* Y1, float * Y2, int W) { 1812 int i = 0; 1813 #ifdef __SSE2__ 1814 const vfloat c1v = F2V(0.2627f); 1815 const vfloat c2v = F2V(0.6780f); 1816 const vfloat c3v = F2V(0.0593f); 1817 for (; i < W - 3; i += 4) { 1818 const vfloat Rv = vmaxf(LVFU(R[i]), ZEROV); 1819 const vfloat Gv = vmaxf(LVFU(G[i]), ZEROV); 1820 const vfloat Bv = vmaxf(LVFU(B[i]), ZEROV); 1821 vfloat yv = c1v * Rv + c2v * Gv + c3v * Bv; 1822 STVFU(Y1[i], yv); 1823 STVFU(Y2[i], yv); 1824 } 1825 #endif 1826 for (; i < W; ++i) { 1827 const float r = std::max(R[i], 0.f); 1828 const float g = std::max(G[i], 0.f); 1829 const float b = std::max(B[i], 0.f); 1830 Y1[i] = Y2[i] = 0.2627f * r + 0.6780f * g + 0.0593f * b; 1831 } 1832 } 1833 1834 }; 1835 1836 } 1837