/src/FreeImage/Source/FreeImage/tmoColorConvert.cpp
C++ | 479 lines | 283 code | 56 blank | 140 comment | 53 complexity | 098073db6c794567be60af307e9df6ad MD5 | raw file
1// ========================================================== 2// High Dynamic Range bitmap conversion routines 3// 4// Design and implementation by 5// - Hervé Drolon (drolon@infonie.fr) 6// - Mihail Naydenov (mnaydenov@users.sourceforge.net) 7// 8// This file is part of FreeImage 3 9// 10// COVERED CODE IS PROVIDED UNDER THIS LICENSE ON AN "AS IS" BASIS, WITHOUT WARRANTY 11// OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES 12// THAT THE COVERED CODE IS FREE OF DEFECTS, MERCHANTABLE, FIT FOR A PARTICULAR PURPOSE 13// OR NON-INFRINGING. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE COVERED 14// CODE IS WITH YOU. SHOULD ANY COVERED CODE PROVE DEFECTIVE IN ANY RESPECT, YOU (NOT 15// THE INITIAL DEVELOPER OR ANY OTHER CONTRIBUTOR) ASSUME THE COST OF ANY NECESSARY 16// SERVICING, REPAIR OR CORRECTION. THIS DISCLAIMER OF WARRANTY CONSTITUTES AN ESSENTIAL 17// PART OF THIS LICENSE. NO USE OF ANY COVERED CODE IS AUTHORIZED HEREUNDER EXCEPT UNDER 18// THIS DISCLAIMER. 19// 20// Use at your own risk! 21// ========================================================== 22 23#include "FreeImage.h" 24#include "Utilities.h" 25#include "ToneMapping.h" 26 27// ---------------------------------------------------------- 28// Convert RGB to and from Yxy, same as in Reinhard et al. SIGGRAPH 2002 29// References : 30// [1] Radiance Home Page [Online] http://radsite.lbl.gov/radiance/HOME.html 31// [2] E. Reinhard, M. Stark, P. Shirley, and J. Ferwerda, 32// Photographic Tone Reproduction for Digital Images, ACM Transactions on Graphics, 33// 21(3):267-276, 2002 (Proceedings of SIGGRAPH 2002). 34// [3] J. Tumblin and H.E. Rushmeier, 35// Tone Reproduction for Realistic Images. IEEE Computer Graphics and Applications, 36// 13(6):42-48, 1993. 37// ---------------------------------------------------------- 38 39/** 40nominal CRT primaries 41*/ 42/* 43static const float CIE_x_r = 0.640F; 44static const float CIE_y_r = 0.330F; 45static const float CIE_x_g = 0.290F; 46static const float CIE_y_g = 0.600F; 47static const float CIE_x_b = 0.150F; 48static const float CIE_y_b = 0.060F; 49static const float CIE_x_w = 0.3333F; // use true white 50static const float CIE_y_w = 0.3333F; 51*/ 52/** 53sRGB primaries 54*/ 55static const float CIE_x_r = 0.640F; 56static const float CIE_y_r = 0.330F; 57static const float CIE_x_g = 0.300F; 58static const float CIE_y_g = 0.600F; 59static const float CIE_x_b = 0.150F; 60static const float CIE_y_b = 0.060F; 61static const float CIE_x_w = 0.3127F; // Illuminant D65 62static const float CIE_y_w = 0.3290F; 63 64static const float CIE_D = ( CIE_x_r*(CIE_y_g - CIE_y_b) + CIE_x_g*(CIE_y_b - CIE_y_r) + CIE_x_b*(CIE_y_r - CIE_y_g) ); 65static const float CIE_C_rD = ( (1/CIE_y_w) * ( CIE_x_w*(CIE_y_g - CIE_y_b) - CIE_y_w*(CIE_x_g - CIE_x_b) + CIE_x_g*CIE_y_b - CIE_x_b*CIE_y_g) ); 66static const float CIE_C_gD = ( (1/CIE_y_w) * ( CIE_x_w*(CIE_y_b - CIE_y_r) - CIE_y_w*(CIE_x_b - CIE_x_r) - CIE_x_r*CIE_y_b + CIE_x_b*CIE_y_r) ); 67static const float CIE_C_bD = ( (1/CIE_y_w) * ( CIE_x_w*(CIE_y_r - CIE_y_g) - CIE_y_w*(CIE_x_r - CIE_x_g) + CIE_x_r*CIE_y_g - CIE_x_g*CIE_y_r) ); 68 69/** 70RGB to XYZ (no white balance) 71*/ 72static const float RGB2XYZ[3][3] = { 73 { CIE_x_r*CIE_C_rD / CIE_D, 74 CIE_x_g*CIE_C_gD / CIE_D, 75 CIE_x_b*CIE_C_bD / CIE_D 76 }, 77 { CIE_y_r*CIE_C_rD / CIE_D, 78 CIE_y_g*CIE_C_gD / CIE_D, 79 CIE_y_b*CIE_C_bD / CIE_D 80 }, 81 { (1 - CIE_x_r-CIE_y_r)*CIE_C_rD / CIE_D, 82 (1 - CIE_x_g-CIE_y_g)*CIE_C_gD / CIE_D, 83 (1 - CIE_x_b-CIE_y_b)*CIE_C_bD / CIE_D 84 } 85}; 86 87/** 88XYZ to RGB (no white balance) 89*/ 90static const float XYZ2RGB[3][3] = { 91 {(CIE_y_g - CIE_y_b - CIE_x_b*CIE_y_g + CIE_y_b*CIE_x_g) / CIE_C_rD, 92 (CIE_x_b - CIE_x_g - CIE_x_b*CIE_y_g + CIE_x_g*CIE_y_b) / CIE_C_rD, 93 (CIE_x_g*CIE_y_b - CIE_x_b*CIE_y_g) / CIE_C_rD 94 }, 95 {(CIE_y_b - CIE_y_r - CIE_y_b*CIE_x_r + CIE_y_r*CIE_x_b) / CIE_C_gD, 96 (CIE_x_r - CIE_x_b - CIE_x_r*CIE_y_b + CIE_x_b*CIE_y_r) / CIE_C_gD, 97 (CIE_x_b*CIE_y_r - CIE_x_r*CIE_y_b) / CIE_C_gD 98 }, 99 {(CIE_y_r - CIE_y_g - CIE_y_r*CIE_x_g + CIE_y_g*CIE_x_r) / CIE_C_bD, 100 (CIE_x_g - CIE_x_r - CIE_x_g*CIE_y_r + CIE_x_r*CIE_y_g) / CIE_C_bD, 101 (CIE_x_r*CIE_y_g - CIE_x_g*CIE_y_r) / CIE_C_bD 102 } 103}; 104 105/** 106This gives approximately the following matrices : 107 108static const float RGB2XYZ[3][3] = { 109 { 0.41239083F, 0.35758433F, 0.18048081F }, 110 { 0.21263903F, 0.71516865F, 0.072192319F }, 111 { 0.019330820F, 0.11919473F, 0.95053220F } 112}; 113static const float XYZ2RGB[3][3] = { 114 { 3.2409699F, -1.5373832F, -0.49861079F }, 115 { -0.96924376F, 1.8759676F, 0.041555084F }, 116 { 0.055630036F, -0.20397687F, 1.0569715F } 117}; 118*/ 119 120// ---------------------------------------------------------- 121 122static const float EPSILON = 1e-06F; 123static const float INF = 1e+10F; 124 125/** 126Convert in-place floating point RGB data to Yxy.<br> 127On output, pixel->red == Y, pixel->green == x, pixel->blue == y 128@param dib Input RGBF / Output Yxy image 129@return Returns TRUE if successful, returns FALSE otherwise 130*/ 131BOOL 132ConvertInPlaceRGBFToYxy(FIBITMAP *dib) { 133 float result[3]; 134 135 if(FreeImage_GetImageType(dib) != FIT_RGBF) 136 return FALSE; 137 138 const unsigned width = FreeImage_GetWidth(dib); 139 const unsigned height = FreeImage_GetHeight(dib); 140 const unsigned pitch = FreeImage_GetPitch(dib); 141 142 BYTE *bits = (BYTE*)FreeImage_GetBits(dib); 143 for(unsigned y = 0; y < height; y++) { 144 FIRGBF *pixel = (FIRGBF*)bits; 145 for(unsigned x = 0; x < width; x++) { 146 result[0] = result[1] = result[2] = 0; 147 for (int i = 0; i < 3; i++) { 148 result[i] += RGB2XYZ[i][0] * pixel[x].red; 149 result[i] += RGB2XYZ[i][1] * pixel[x].green; 150 result[i] += RGB2XYZ[i][2] * pixel[x].blue; 151 } 152 const float W = result[0] + result[1] + result[2]; 153 const float Y = result[1]; 154 if(W > 0) { 155 pixel[x].red = Y; // Y 156 pixel[x].green = result[0] / W; // x 157 pixel[x].blue = result[1] / W; // y 158 } else { 159 pixel[x].red = pixel[x].green = pixel[x].blue = 0; 160 } 161 } 162 // next line 163 bits += pitch; 164 } 165 166 return TRUE; 167} 168 169/** 170Convert in-place Yxy image to floating point RGB data.<br> 171On input, pixel->red == Y, pixel->green == x, pixel->blue == y 172@param dib Input Yxy / Output RGBF image 173@return Returns TRUE if successful, returns FALSE otherwise 174*/ 175BOOL 176ConvertInPlaceYxyToRGBF(FIBITMAP *dib) { 177 float result[3]; 178 float X, Y, Z; 179 180 if(FreeImage_GetImageType(dib) != FIT_RGBF) 181 return FALSE; 182 183 const unsigned width = FreeImage_GetWidth(dib); 184 const unsigned height = FreeImage_GetHeight(dib); 185 const unsigned pitch = FreeImage_GetPitch(dib); 186 187 BYTE *bits = (BYTE*)FreeImage_GetBits(dib); 188 for(unsigned y = 0; y < height; y++) { 189 FIRGBF *pixel = (FIRGBF*)bits; 190 for(unsigned x = 0; x < width; x++) { 191 Y = pixel[x].red; // Y 192 result[1] = pixel[x].green; // x 193 result[2] = pixel[x].blue; // y 194 if ((Y > EPSILON) && (result[1] > EPSILON) && (result[2] > EPSILON)) { 195 X = (result[1] * Y) / result[2]; 196 Z = (X / result[1]) - X - Y; 197 } else { 198 X = Z = EPSILON; 199 } 200 pixel[x].red = X; 201 pixel[x].green = Y; 202 pixel[x].blue = Z; 203 result[0] = result[1] = result[2] = 0; 204 for (int i = 0; i < 3; i++) { 205 result[i] += XYZ2RGB[i][0] * pixel[x].red; 206 result[i] += XYZ2RGB[i][1] * pixel[x].green; 207 result[i] += XYZ2RGB[i][2] * pixel[x].blue; 208 } 209 pixel[x].red = result[0]; // R 210 pixel[x].green = result[1]; // G 211 pixel[x].blue = result[2]; // B 212 } 213 // next line 214 bits += pitch; 215 } 216 217 return TRUE; 218} 219 220/** 221Get the maximum, minimum and average luminance.<br> 222On input, pixel->red == Y, pixel->green == x, pixel->blue == y 223@param Yxy Source Yxy image to analyze 224@param maxLum Maximum luminance 225@param minLum Minimum luminance 226@param worldLum Average luminance (world adaptation luminance) 227@return Returns TRUE if successful, returns FALSE otherwise 228*/ 229BOOL 230LuminanceFromYxy(FIBITMAP *Yxy, float *maxLum, float *minLum, float *worldLum) { 231 if(FreeImage_GetImageType(Yxy) != FIT_RGBF) 232 return FALSE; 233 234 const unsigned width = FreeImage_GetWidth(Yxy); 235 const unsigned height = FreeImage_GetHeight(Yxy); 236 const unsigned pitch = FreeImage_GetPitch(Yxy); 237 238 float max_lum = 0, min_lum = 0; 239 double sum = 0; 240 241 BYTE *bits = (BYTE*)FreeImage_GetBits(Yxy); 242 for(unsigned y = 0; y < height; y++) { 243 const FIRGBF *pixel = (FIRGBF*)bits; 244 for(unsigned x = 0; x < width; x++) { 245 const float Y = pixel[x].red; 246 max_lum = (max_lum < Y) ? Y : max_lum; // max Luminance in the scene 247 min_lum = (min_lum < Y) ? min_lum : Y; // min Luminance in the scene 248 sum += log(2.3e-5F + Y); // contrast constant in Tumblin paper 249 } 250 // next line 251 bits += pitch; 252 } 253 // maximum luminance 254 *maxLum = max_lum; 255 // minimum luminance 256 *minLum = min_lum; 257 // average log luminance 258 double avgLogLum = (sum / (width * height)); 259 // world adaptation luminance 260 *worldLum = (float)exp(avgLogLum); 261 262 return TRUE; 263} 264 265/** 266Clamp RGBF image highest values to display white, 267then convert to 24-bit RGB 268*/ 269FIBITMAP* 270ClampConvertRGBFTo24(FIBITMAP *src) { 271 if(FreeImage_GetImageType(src) != FIT_RGBF) 272 return FALSE; 273 274 const unsigned width = FreeImage_GetWidth(src); 275 const unsigned height = FreeImage_GetHeight(src); 276 277 FIBITMAP *dst = FreeImage_Allocate(width, height, 24, FI_RGBA_RED_MASK, FI_RGBA_GREEN_MASK, FI_RGBA_BLUE_MASK); 278 if(!dst) return NULL; 279 280 const unsigned src_pitch = FreeImage_GetPitch(src); 281 const unsigned dst_pitch = FreeImage_GetPitch(dst); 282 283 BYTE *src_bits = (BYTE*)FreeImage_GetBits(src); 284 BYTE *dst_bits = (BYTE*)FreeImage_GetBits(dst); 285 286 for(unsigned y = 0; y < height; y++) { 287 const FIRGBF *src_pixel = (FIRGBF*)src_bits; 288 BYTE *dst_pixel = (BYTE*)dst_bits; 289 for(unsigned x = 0; x < width; x++) { 290 const float red = (src_pixel[x].red > 1) ? 1 : src_pixel[x].red; 291 const float green = (src_pixel[x].green > 1) ? 1 : src_pixel[x].green; 292 const float blue = (src_pixel[x].blue > 1) ? 1 : src_pixel[x].blue; 293 294 dst_pixel[FI_RGBA_RED] = (BYTE)(255.0F * red + 0.5F); 295 dst_pixel[FI_RGBA_GREEN] = (BYTE)(255.0F * green + 0.5F); 296 dst_pixel[FI_RGBA_BLUE] = (BYTE)(255.0F * blue + 0.5F); 297 dst_pixel += 3; 298 } 299 src_bits += src_pitch; 300 dst_bits += dst_pitch; 301 } 302 303 return dst; 304} 305 306/** 307Extract the luminance channel L from a RGBF image. 308Luminance is calculated from the sRGB model (RGB2XYZ matrix) 309using a D65 white point : 310L = ( 0.2126 * r ) + ( 0.7152 * g ) + ( 0.0722 * b ) 311Reference : 312A Standard Default Color Space for the Internet - sRGB. 313[online] http://www.w3.org/Graphics/Color/sRGB 314*/ 315FIBITMAP* 316ConvertRGBFToY(FIBITMAP *src) { 317 if(FreeImage_GetImageType(src) != FIT_RGBF) 318 return FALSE; 319 320 const unsigned width = FreeImage_GetWidth(src); 321 const unsigned height = FreeImage_GetHeight(src); 322 323 FIBITMAP *dst = FreeImage_AllocateT(FIT_FLOAT, width, height); 324 if(!dst) return NULL; 325 326 const unsigned src_pitch = FreeImage_GetPitch(src); 327 const unsigned dst_pitch = FreeImage_GetPitch(dst); 328 329 330 BYTE *src_bits = (BYTE*)FreeImage_GetBits(src); 331 BYTE *dst_bits = (BYTE*)FreeImage_GetBits(dst); 332 333 for(unsigned y = 0; y < height; y++) { 334 const FIRGBF *src_pixel = (FIRGBF*)src_bits; 335 float *dst_pixel = (float*)dst_bits; 336 for(unsigned x = 0; x < width; x++) { 337 const float L = LUMA_REC709(src_pixel[x].red, src_pixel[x].green, src_pixel[x].blue); 338 dst_pixel[x] = (L > 0) ? L : 0; 339 } 340 // next line 341 src_bits += src_pitch; 342 dst_bits += dst_pitch; 343 } 344 345 return dst; 346} 347 348/** 349Get the maximum, minimum, average luminance and log average luminance from a Y image 350@param dib Source Y image to analyze 351@param maxLum Maximum luminance 352@param minLum Minimum luminance 353@param Lav Average luminance 354@param Llav Log average luminance (also known as 'world adaptation luminance') 355@return Returns TRUE if successful, returns FALSE otherwise 356@see ConvertRGBFToY, FreeImage_TmoReinhard05Ex 357*/ 358BOOL 359LuminanceFromY(FIBITMAP *dib, float *maxLum, float *minLum, float *Lav, float *Llav) { 360 if(FreeImage_GetImageType(dib) != FIT_FLOAT) 361 return FALSE; 362 363 unsigned width = FreeImage_GetWidth(dib); 364 unsigned height = FreeImage_GetHeight(dib); 365 unsigned pitch = FreeImage_GetPitch(dib); 366 367 float max_lum = -1e20F, min_lum = 1e20F; 368 double sumLum = 0, sumLogLum = 0; 369 370 BYTE *bits = (BYTE*)FreeImage_GetBits(dib); 371 for(unsigned y = 0; y < height; y++) { 372 const float *pixel = (float*)bits; 373 for(unsigned x = 0; x < width; x++) { 374 const float Y = pixel[x]; 375 max_lum = (max_lum < Y) ? Y : max_lum; // max Luminance in the scene 376 min_lum = ((Y > 0) && (min_lum < Y)) ? min_lum : Y; // min Luminance in the scene 377 sumLum += Y; // average luminance 378 sumLogLum += log(2.3e-5F + Y); // contrast constant in Tumblin paper 379 } 380 // next line 381 bits += pitch; 382 } 383 384 // maximum luminance 385 *maxLum = max_lum; 386 // minimum luminance 387 *minLum = min_lum; 388 // average luminance 389 *Lav = (float)(sumLum / (width * height)); 390 // average log luminance, a.k.a. world adaptation luminance 391 *Llav = (float)exp(sumLogLum / (width * height)); 392 393 return TRUE; 394} 395// -------------------------------------------------------------------------- 396 397static void findMaxMinPercentile(FIBITMAP *Y, float minPrct, float *minLum, float maxPrct, float *maxLum) { 398 int x, y; 399 int width = FreeImage_GetWidth(Y); 400 int height = FreeImage_GetHeight(Y); 401 int pitch = FreeImage_GetPitch(Y); 402 403 std::vector<float> vY(width * height); 404 405 BYTE *bits = (BYTE*)FreeImage_GetBits(Y); 406 for(y = 0; y < height; y++) { 407 float *pixel = (float*)bits; 408 for(x = 0; x < width; x++) { 409 if(pixel[x] != 0) { 410 vY.push_back(pixel[x]); 411 } 412 } 413 bits += pitch; 414 } 415 416 std::sort(vY.begin(), vY.end()); 417 418 *minLum = vY.at( int(minPrct * vY.size()) ); 419 *maxLum = vY.at( int(maxPrct * vY.size()) ); 420} 421 422/** 423Clipping function<br> 424Remove any extremely bright and/or extremely dark pixels 425and normalize between 0 and 1. 426@param Y Input/Output image 427@param minPrct Minimum percentile 428@param maxPrct Maximum percentile 429*/ 430void 431NormalizeY(FIBITMAP *Y, float minPrct, float maxPrct) { 432 int x, y; 433 float maxLum, minLum; 434 435 if(minPrct > maxPrct) { 436 // swap values 437 float t = minPrct; minPrct = maxPrct; maxPrct = t; 438 } 439 if(minPrct < 0) minPrct = 0; 440 if(maxPrct > 1) maxPrct = 1; 441 442 int width = FreeImage_GetWidth(Y); 443 int height = FreeImage_GetHeight(Y); 444 int pitch = FreeImage_GetPitch(Y); 445 446 // find max & min luminance values 447 if((minPrct > 0) || (maxPrct < 1)) { 448 maxLum = 0, minLum = 0; 449 findMaxMinPercentile(Y, minPrct, &minLum, maxPrct, &maxLum); 450 } else { 451 maxLum = -1e20F, minLum = 1e20F; 452 BYTE *bits = (BYTE*)FreeImage_GetBits(Y); 453 for(y = 0; y < height; y++) { 454 const float *pixel = (float*)bits; 455 for(x = 0; x < width; x++) { 456 const float value = pixel[x]; 457 maxLum = (maxLum < value) ? value : maxLum; // max Luminance in the scene 458 minLum = (minLum < value) ? minLum : value; // min Luminance in the scene 459 } 460 // next line 461 bits += pitch; 462 } 463 } 464 if(maxLum == minLum) return; 465 466 // normalize to range 0..1 467 const float divider = maxLum - minLum; 468 BYTE *bits = (BYTE*)FreeImage_GetBits(Y); 469 for(y = 0; y < height; y++) { 470 float *pixel = (float*)bits; 471 for(x = 0; x < width; x++) { 472 pixel[x] = (pixel[x] - minLum) / divider; 473 if(pixel[x] <= 0) pixel[x] = EPSILON; 474 if(pixel[x] > 1) pixel[x] = 1; 475 } 476 // next line 477 bits += pitch; 478 } 479}