/src/FreeImage/Source/FreeImage/tmoColorConvert.cpp

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