/src/FreeImage/Source/FreeImage/tmoDrago03.cpp

https://bitbucket.org/cabalistic/ogredeps/ · C++ · 295 lines · 166 code · 38 blank · 91 comment · 36 complexity · 122d9ef476c636a6137b8b9fcc664ada MD5 · raw file

  1. // ==========================================================
  2. // Tone mapping operator (Drago, 2003)
  3. //
  4. // Design and implementation by
  5. // - Hervé Drolon (drolon@infonie.fr)
  6. //
  7. // This file is part of FreeImage 3
  8. //
  9. // COVERED CODE IS PROVIDED UNDER THIS LICENSE ON AN "AS IS" BASIS, WITHOUT WARRANTY
  10. // OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES
  11. // THAT THE COVERED CODE IS FREE OF DEFECTS, MERCHANTABLE, FIT FOR A PARTICULAR PURPOSE
  12. // OR NON-INFRINGING. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE COVERED
  13. // CODE IS WITH YOU. SHOULD ANY COVERED CODE PROVE DEFECTIVE IN ANY RESPECT, YOU (NOT
  14. // THE INITIAL DEVELOPER OR ANY OTHER CONTRIBUTOR) ASSUME THE COST OF ANY NECESSARY
  15. // SERVICING, REPAIR OR CORRECTION. THIS DISCLAIMER OF WARRANTY CONSTITUTES AN ESSENTIAL
  16. // PART OF THIS LICENSE. NO USE OF ANY COVERED CODE IS AUTHORIZED HEREUNDER EXCEPT UNDER
  17. // THIS DISCLAIMER.
  18. //
  19. // Use at your own risk!
  20. // ==========================================================
  21. #include "FreeImage.h"
  22. #include "Utilities.h"
  23. #include "ToneMapping.h"
  24. // ----------------------------------------------------------
  25. // Logarithmic mapping operator
  26. // Reference:
  27. // [1] F. Drago, K. Myszkowski, T. Annen, and N. Chiba,
  28. // Adaptive Logarithmic Mapping for Displaying High Contrast Scenes,
  29. // Eurographics 2003.
  30. // ----------------------------------------------------------
  31. /**
  32. Bias function
  33. */
  34. static inline double
  35. biasFunction(const double b, const double x) {
  36. return pow (x, b); // pow(x, log(bias)/log(0.5)
  37. }
  38. /**
  39. Padé approximation of log(x + 1)
  40. x(6+x)/(6+4x) good if x < 1
  41. x*(6 + 0.7662x)/(5.9897 + 3.7658x) between 1 and 2
  42. See http://www.nezumi.demon.co.uk/consult/logx.htm
  43. */
  44. static inline double
  45. pade_log(const double x) {
  46. if(x < 1) {
  47. return (x * (6 + x) / (6 + 4 * x));
  48. } else if(x < 2) {
  49. return (x * (6 + 0.7662 * x) / (5.9897 + 3.7658 * x));
  50. }
  51. return log(x + 1);
  52. }
  53. /**
  54. Log mapping operator
  55. @param dib Input / Output Yxy image
  56. @param maxLum Maximum luminance
  57. @param avgLum Average luminance (world adaptation luminance)
  58. @param biasParam Bias parameter (a zero value default to 0.85)
  59. @param exposure Exposure parameter (default to 0)
  60. @return Returns TRUE if successful, returns FALSE otherwise
  61. */
  62. static BOOL
  63. ToneMappingDrago03(FIBITMAP *dib, const float maxLum, const float avgLum, float biasParam, const float exposure) {
  64. const float LOG05 = -0.693147F; // log(0.5)
  65. double Lmax, divider, interpol, biasP;
  66. unsigned x, y;
  67. double L;
  68. if(FreeImage_GetImageType(dib) != FIT_RGBF)
  69. return FALSE;
  70. const unsigned width = FreeImage_GetWidth(dib);
  71. const unsigned height = FreeImage_GetHeight(dib);
  72. const unsigned pitch = FreeImage_GetPitch(dib);
  73. // arbitrary Bias Parameter
  74. if(biasParam == 0)
  75. biasParam = 0.85F;
  76. // normalize maximum luminance by average luminance
  77. Lmax = maxLum / avgLum;
  78. divider = log10(Lmax+1);
  79. biasP = log(biasParam)/LOG05;
  80. #if !defined(DRAGO03_FAST)
  81. /**
  82. Normal tone mapping of every pixel
  83. further acceleration is obtained by a Padé approximation of log(x + 1)
  84. */
  85. BYTE *bits = (BYTE*)FreeImage_GetBits(dib);
  86. for(y = 0; y < height; y++) {
  87. FIRGBF *pixel = (FIRGBF*)bits;
  88. for(x = 0; x < width; x++) {
  89. double Yw = pixel[x].red / avgLum;
  90. Yw *= exposure;
  91. interpol = log(2 + biasFunction(biasP, Yw / Lmax) * 8);
  92. L = pade_log(Yw);// log(Yw + 1)
  93. pixel[x].red = (float)((L / interpol) / divider);
  94. }
  95. // next line
  96. bits += pitch;
  97. }
  98. #else
  99. unsigned index;
  100. int i, j;
  101. unsigned max_width = width - (width % 3);
  102. unsigned max_height = height - (height % 3);
  103. unsigned fpitch = pitch / sizeof(FIRGBF);
  104. /**
  105. fast tone mapping
  106. split the image into 3x3 pixel tiles and perform the computation for each group of 9 pixels
  107. further acceleration is obtained by a Padé approximation of log(x + 1)
  108. => produce artifacts and not so faster, so the code has been disabled
  109. */
  110. #define PIXEL(x, y) image[y*fpitch + x].red
  111. FIRGBF *image = (FIRGBF*)FreeImage_GetBits(dib);
  112. for(y = 0; y < max_height; y += 3) {
  113. for(x = 0; x < max_width; x += 3) {
  114. double average = 0;
  115. for(i = 0; i < 3; i++) {
  116. for(j = 0; j < 3; j++) {
  117. index = (y + i)*fpitch + (x + j);
  118. image[index].red /= (float)avgLum;
  119. image[index].red *= exposure;
  120. average += image[index].red;
  121. }
  122. }
  123. average = average / 9 - PIXEL(x, y);
  124. if(average > -1 && average < 1) {
  125. interpol = log(2 + pow(PIXEL(x + 1, y + 1) / Lmax, biasP) * 8);
  126. for(i = 0; i < 3; i++) {
  127. for(j = 0; j < 3; j++) {
  128. index = (y + i)*fpitch + (x + j);
  129. L = pade_log(image[index].red);// log(image[index].red + 1)
  130. image[index].red = (float)((L / interpol) / divider);
  131. }
  132. }
  133. }
  134. else {
  135. for(i = 0; i < 3; i++) {
  136. for(j = 0; j < 3; j++) {
  137. index = (y + i)*fpitch + (x + j);
  138. interpol = log(2 + pow(image[index].red / Lmax, biasP) * 8);
  139. L = pade_log(image[index].red);// log(image[index].red + 1)
  140. image[index].red = (float)((L / interpol) / divider);
  141. }
  142. }
  143. }
  144. } //x
  145. } // y
  146. /**
  147. Normal tone mapping of every pixel for the remaining right and bottom bands
  148. */
  149. BYTE *bits;
  150. // right band
  151. bits = (BYTE*)FreeImage_GetBits(dib);
  152. for(y = 0; y < height; y++) {
  153. FIRGBF *pixel = (FIRGBF*)bits;
  154. for(x = max_width; x < width; x++) {
  155. double Yw = pixel[x].red / avgLum;
  156. Yw *= exposure;
  157. interpol = log(2 + biasFunction(biasP, Yw / Lmax) * 8);
  158. L = pade_log(Yw);// log(Yw + 1)
  159. pixel[x].red = (float)((L / interpol) / divider);
  160. }
  161. // next line
  162. bits += pitch;
  163. }
  164. // bottom band
  165. bits = (BYTE*)FreeImage_GetBits(dib);
  166. for(y = max_height; y < height; y++) {
  167. FIRGBF *pixel = (FIRGBF*)bits;
  168. for(x = 0; x < max_width; x++) {
  169. double Yw = pixel[x].red / avgLum;
  170. Yw *= exposure;
  171. interpol = log(2 + biasFunction(biasP, Yw / Lmax) * 8);
  172. L = pade_log(Yw);// log(Yw + 1)
  173. pixel[x].red = (float)((L / interpol) / divider);
  174. }
  175. // next line
  176. bits += pitch;
  177. }
  178. #endif // DRAGO03_FAST
  179. return TRUE;
  180. }
  181. /**
  182. Custom gamma correction based on the ITU-R BT.709 standard
  183. @param dib RGBF image to be corrected
  184. @param gammaval Gamma value (2.2 is a good default value)
  185. @return Returns TRUE if successful, returns FALSE otherwise
  186. */
  187. static BOOL
  188. REC709GammaCorrection(FIBITMAP *dib, const float gammaval) {
  189. if(FreeImage_GetImageType(dib) != FIT_RGBF)
  190. return FALSE;
  191. float slope = 4.5F;
  192. float start = 0.018F;
  193. const float fgamma = (float)((0.45 / gammaval) * 2);
  194. if(gammaval >= 2.1F) {
  195. start = (float)(0.018 / ((gammaval - 2) * 7.5));
  196. slope = (float)(4.5 * ((gammaval - 2) * 7.5));
  197. } else if (gammaval <= 1.9F) {
  198. start = (float)(0.018 * ((2 - gammaval) * 7.5));
  199. slope = (float)(4.5 / ((2 - gammaval) * 7.5));
  200. }
  201. const unsigned width = FreeImage_GetWidth(dib);
  202. const unsigned height = FreeImage_GetHeight(dib);
  203. const unsigned pitch = FreeImage_GetPitch(dib);
  204. BYTE *bits = (BYTE*)FreeImage_GetBits(dib);
  205. for(unsigned y = 0; y < height; y++) {
  206. float *pixel = (float*)bits;
  207. for(unsigned x = 0; x < width; x++) {
  208. for(int i = 0; i < 3; i++) {
  209. *pixel = (*pixel <= start) ? *pixel * slope : (1.099F * pow(*pixel, fgamma) - 0.099F);
  210. pixel++;
  211. }
  212. }
  213. bits += pitch;
  214. }
  215. return TRUE;
  216. }
  217. // ----------------------------------------------------------
  218. // Main algorithm
  219. // ----------------------------------------------------------
  220. /**
  221. Apply the Adaptive Logarithmic Mapping operator to a HDR image and convert to 24-bit RGB
  222. @param src Input RGB16 or RGB[A]F image
  223. @param gamma Gamma correction (gamma > 0). 1 means no correction, 2.2 in the original paper.
  224. @param exposure Exposure parameter (0 means no correction, 0 in the original paper)
  225. @return Returns a 24-bit RGB image if successful, returns NULL otherwise
  226. */
  227. FIBITMAP* DLL_CALLCONV
  228. FreeImage_TmoDrago03(FIBITMAP *src, double gamma, double exposure) {
  229. float maxLum, minLum, avgLum;
  230. if(!FreeImage_HasPixels(src)) return NULL;
  231. // working RGBF variable
  232. FIBITMAP *dib = NULL;
  233. dib = FreeImage_ConvertToRGBF(src);
  234. if(!dib) return NULL;
  235. // default algorithm parameters
  236. const float biasParam = 0.85F;
  237. const float expoParam = (float)pow(2.0, exposure); //default exposure is 1, 2^0
  238. // convert to Yxy
  239. ConvertInPlaceRGBFToYxy(dib);
  240. // get the luminance
  241. LuminanceFromYxy(dib, &maxLum, &minLum, &avgLum);
  242. // perform the tone mapping
  243. ToneMappingDrago03(dib, maxLum, avgLum, biasParam, expoParam);
  244. // convert back to RGBF
  245. ConvertInPlaceYxyToRGBF(dib);
  246. if(gamma != 1) {
  247. // perform gamma correction
  248. REC709GammaCorrection(dib, (float)gamma);
  249. }
  250. // clamp image highest values to display white, then convert to 24-bit RGB
  251. FIBITMAP *dst = ClampConvertRGBFTo24(dib);
  252. // clean-up and return
  253. FreeImage_Unload(dib);
  254. // copy metadata from src to dst
  255. FreeImage_CloneMetadata(dst, src);
  256. return dst;
  257. }