1 // ==========================================================
2 // Tone mapping operator (Fattal, 2002)
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
22 #include "FreeImage.h"
23 #include "Utilities.h"
24 #include "ToneMapping.h"
25
26 // ----------------------------------------------------------
27 // Gradient domain HDR compression
28 // Reference:
29 // [1] R. Fattal, D. Lischinski, and M.Werman,
30 // Gradient domain high dynamic range compression,
31 // ACM Transactions on Graphics, special issue on Proc. of ACM SIGGRAPH 2002,
32 // San Antonio, Texas, vol. 21(3), pp. 257-266, 2002.
33 // ----------------------------------------------------------
34
35 static const float EPSILON = 1e-4F;
36
37 /**
38 Performs a 5 by 5 gaussian filtering using two 1D convolutions,
39 followed by a subsampling by 2.
40 @param dib Input image
41 @return Returns a blurred image of size SIZE(dib)/2
42 @see GaussianPyramid
43 */
GaussianLevel5x5(FIBITMAP * dib)44 static FIBITMAP* GaussianLevel5x5(FIBITMAP *dib) {
45 FIBITMAP *h_dib = NULL, *v_dib = NULL, *dst = NULL;
46 float *src_pixel, *dst_pixel;
47
48 try {
49 const FREE_IMAGE_TYPE image_type = FreeImage_GetImageType(dib);
50 if(image_type != FIT_FLOAT) throw(1);
51
52 const unsigned width = FreeImage_GetWidth(dib);
53 const unsigned height = FreeImage_GetHeight(dib);
54
55 h_dib = FreeImage_AllocateT(image_type, width, height);
56 v_dib = FreeImage_AllocateT(image_type, width, height);
57 if(!h_dib || !v_dib) throw(1);
58
59 const unsigned pitch = FreeImage_GetPitch(dib) / sizeof(float);
60
61 // horizontal convolution dib -> h_dib
62
63 src_pixel = (float*)FreeImage_GetBits(dib);
64 dst_pixel = (float*)FreeImage_GetBits(h_dib);
65
66 for(unsigned y = 0; y < height; y++) {
67 // work on line y
68 for(unsigned x = 2; x < width - 2; x++) {
69 dst_pixel[x] = src_pixel[x-2] + src_pixel[x+2] + 4 * (src_pixel[x-1] + src_pixel[x+1]) + 6 * src_pixel[x];
70 dst_pixel[x] /= 16;
71 }
72 // boundary mirroring
73 dst_pixel[0] = (2 * src_pixel[2] + 8 * src_pixel[1] + 6 * src_pixel[0]) / 16;
74 dst_pixel[1] = (src_pixel[3] + 4 * (src_pixel[0] + src_pixel[2]) + 7 * src_pixel[1]) / 16;
75 dst_pixel[width-2] = (src_pixel[width-4] + 5 * src_pixel[width-1] + 4 * src_pixel[width-3] + 6 * src_pixel[width-2]) / 16;
76 dst_pixel[width-1] = (src_pixel[width-3] + 5 * src_pixel[width-2] + 10 * src_pixel[width-1]) / 16;
77
78 // next line
79 src_pixel += pitch;
80 dst_pixel += pitch;
81 }
82
83 // vertical convolution h_dib -> v_dib
84
85 src_pixel = (float*)FreeImage_GetBits(h_dib);
86 dst_pixel = (float*)FreeImage_GetBits(v_dib);
87
88 for(unsigned x = 0; x < width; x++) {
89 // work on column x
90 for(unsigned y = 2; y < height - 2; y++) {
91 const unsigned index = y*pitch + x;
92 dst_pixel[index] = src_pixel[index-2*pitch] + src_pixel[index+2*pitch] + 4 * (src_pixel[index-pitch] + src_pixel[index+pitch]) + 6 * src_pixel[index];
93 dst_pixel[index] /= 16;
94 }
95 // boundary mirroring
96 dst_pixel[x] = (2 * src_pixel[x+2*pitch] + 8 * src_pixel[x+pitch] + 6 * src_pixel[x]) / 16;
97 dst_pixel[x+pitch] = (src_pixel[x+3*pitch] + 4 * (src_pixel[x] + src_pixel[x+2*pitch]) + 7 * src_pixel[x+pitch]) / 16;
98 dst_pixel[(height-2)*pitch+x] = (src_pixel[(height-4)*pitch+x] + 5 * src_pixel[(height-1)*pitch+x] + 4 * src_pixel[(height-3)*pitch+x] + 6 * src_pixel[(height-2)*pitch+x]) / 16;
99 dst_pixel[(height-1)*pitch+x] = (src_pixel[(height-3)*pitch+x] + 5 * src_pixel[(height-2)*pitch+x] + 10 * src_pixel[(height-1)*pitch+x]) / 16;
100 }
101
102 FreeImage_Unload(h_dib); h_dib = NULL;
103
104 // perform downsampling
105
106 dst = FreeImage_Rescale(v_dib, width/2, height/2, FILTER_BILINEAR);
107
108 FreeImage_Unload(v_dib);
109
110 return dst;
111
112 } catch(int) {
113 if(h_dib) FreeImage_Unload(h_dib);
114 if(v_dib) FreeImage_Unload(v_dib);
115 if(dst) FreeImage_Unload(dst);
116 return NULL;
117 }
118 }
119
120 /**
121 Compute a Gaussian pyramid using the specified number of levels.
122 @param H Original bitmap
123 @param pyramid Resulting pyramid array
124 @param nlevels Number of resolution levels
125 @return Returns TRUE if successful, returns FALSE otherwise
126 */
GaussianPyramid(FIBITMAP * H,FIBITMAP ** pyramid,int nlevels)127 static BOOL GaussianPyramid(FIBITMAP *H, FIBITMAP **pyramid, int nlevels) {
128 try {
129 // first level is the original image
130 pyramid[0] = FreeImage_Clone(H);
131 if(pyramid[0] == NULL) throw(1);
132 // compute next levels
133 for(int k = 1; k < nlevels; k++) {
134 pyramid[k] = GaussianLevel5x5(pyramid[k-1]);
135 if(pyramid[k] == NULL) throw(1);
136 }
137 return TRUE;
138 } catch(int) {
139 for(int k = 0; k < nlevels; k++) {
140 if(pyramid[k] != NULL) {
141 FreeImage_Unload(pyramid[k]);
142 pyramid[k] = NULL;
143 }
144 }
145 return FALSE;
146 }
147 }
148
149 /**
150 Compute the gradient magnitude of an input image H using central differences,
151 and returns the average gradient.
152 @param H Input image
153 @param avgGrad [out] Average gradient
154 @param k Level number
155 @return Returns the gradient magnitude if successful, returns NULL otherwise
156 @see GradientPyramid
157 */
GradientLevel(FIBITMAP * H,float * avgGrad,int k)158 static FIBITMAP* GradientLevel(FIBITMAP *H, float *avgGrad, int k) {
159 FIBITMAP *G = NULL;
160
161 try {
162 const FREE_IMAGE_TYPE image_type = FreeImage_GetImageType(H);
163 if(image_type != FIT_FLOAT) throw(1);
164
165 const unsigned width = FreeImage_GetWidth(H);
166 const unsigned height = FreeImage_GetHeight(H);
167
168 G = FreeImage_AllocateT(image_type, width, height);
169 if(!G) throw(1);
170
171 const unsigned pitch = FreeImage_GetPitch(H) / sizeof(float);
172
173 const float divider = (float)(1 << (k + 1));
174 float average = 0;
175
176 float *src_pixel = (float*)FreeImage_GetBits(H);
177 float *dst_pixel = (float*)FreeImage_GetBits(G);
178
179 for(unsigned y = 0; y < height; y++) {
180 const unsigned n = (y == 0 ? 0 : y-1);
181 const unsigned s = (y+1 == height ? y : y+1);
182 for(unsigned x = 0; x < width; x++) {
183 const unsigned w = (x == 0 ? 0 : x-1);
184 const unsigned e = (x+1 == width ? x : x+1);
185 // central difference
186 const float gx = (src_pixel[y*pitch+e] - src_pixel[y*pitch+w]) / divider; // [Hk(x+1, y) - Hk(x-1, y)] / 2**(k+1)
187 const float gy = (src_pixel[s*pitch+x] - src_pixel[n*pitch+x]) / divider; // [Hk(x, y+1) - Hk(x, y-1)] / 2**(k+1)
188 // gradient
189 dst_pixel[x] = sqrt(gx*gx + gy*gy);
190 // average gradient
191 average += dst_pixel[x];
192 }
193 // next line
194 dst_pixel += pitch;
195 }
196
197 *avgGrad = average / (width * height);
198
199 return G;
200
201 } catch(int) {
202 if(G) FreeImage_Unload(G);
203 return NULL;
204 }
205 }
206
207 /**
208 Calculate gradient magnitude and its average value on each pyramid level
209 @param pyramid Gaussian pyramid (nlevels levels)
210 @param nlevels Number of levels
211 @param gradients [out] Gradient pyramid (nlevels levels)
212 @param avgGrad [out] Average gradient on each level (array of size nlevels)
213 @return Returns TRUE if successful, returns FALSE otherwise
214 */
GradientPyramid(FIBITMAP ** pyramid,int nlevels,FIBITMAP ** gradients,float * avgGrad)215 static BOOL GradientPyramid(FIBITMAP **pyramid, int nlevels, FIBITMAP **gradients, float *avgGrad) {
216 try {
217 for(int k = 0; k < nlevels; k++) {
218 FIBITMAP *Hk = pyramid[k];
219 gradients[k] = GradientLevel(Hk, &avgGrad[k], k);
220 if(gradients[k] == NULL) throw(1);
221 }
222 return TRUE;
223 } catch(int) {
224 for(int k = 0; k < nlevels; k++) {
225 if(gradients[k] != NULL) {
226 FreeImage_Unload(gradients[k]);
227 gradients[k] = NULL;
228 }
229 }
230 return FALSE;
231 }
232 }
233
234 /**
235 Compute the gradient attenuation function PHI(x, y)
236 @param gradients Gradient pyramid (nlevels levels)
237 @param avgGrad Average gradient on each level (array of size nlevels)
238 @param nlevels Number of levels
239 @param alpha Parameter alpha in the paper
240 @param beta Parameter beta in the paper
241 @return Returns the attenuation matrix Phi if successful, returns NULL otherwise
242 */
PhiMatrix(FIBITMAP ** gradients,float * avgGrad,int nlevels,float alpha,float beta)243 static FIBITMAP* PhiMatrix(FIBITMAP **gradients, float *avgGrad, int nlevels, float alpha, float beta) {
244 float *src_pixel, *dst_pixel;
245 FIBITMAP **phi = NULL;
246
247 try {
248 phi = (FIBITMAP**)malloc(nlevels * sizeof(FIBITMAP*));
249 if(!phi) throw(1);
250 memset(phi, 0, nlevels * sizeof(FIBITMAP*));
251
252 for(int k = nlevels-1; k >= 0; k--) {
253 // compute phi(k)
254
255 FIBITMAP *Gk = gradients[k];
256
257 const unsigned width = FreeImage_GetWidth(Gk);
258 const unsigned height = FreeImage_GetHeight(Gk);
259 const unsigned pitch = FreeImage_GetPitch(Gk) / sizeof(float);
260
261 // parameter alpha is 0.1 times the average gradient magnitude
262 // also, note the factor of 2**k in the denominator;
263 // that is there to correct for the fact that an average gradient avgGrad(H) over 2**k pixels
264 // in the original image will appear as a gradient grad(Hk) = 2**k*avgGrad(H) over a single pixel in Hk.
265 float ALPHA = alpha * avgGrad[k] * (float)((int)1 << k);
266 if(ALPHA == 0) ALPHA = EPSILON;
267
268 phi[k] = FreeImage_AllocateT(FIT_FLOAT, width, height);
269 if(!phi[k]) throw(1);
270
271 src_pixel = (float*)FreeImage_GetBits(Gk);
272 dst_pixel = (float*)FreeImage_GetBits(phi[k]);
273 for(unsigned y = 0; y < height; y++) {
274 for(unsigned x = 0; x < width; x++) {
275 // compute (alpha / grad) * (grad / alpha) ** beta
276 const float v = src_pixel[x] / ALPHA;
277 const float value = (float)pow((float)v, (float)(beta-1));
278 dst_pixel[x] = (value > 1) ? 1 : value;
279 }
280 // next line
281 src_pixel += pitch;
282 dst_pixel += pitch;
283 }
284
285 if(k < nlevels-1) {
286 // compute PHI(k) = L( PHI(k+1) ) * phi(k)
287 FIBITMAP *L = FreeImage_Rescale(phi[k+1], width, height, FILTER_BILINEAR);
288 if(!L) throw(1);
289
290 src_pixel = (float*)FreeImage_GetBits(L);
291 dst_pixel = (float*)FreeImage_GetBits(phi[k]);
292 for(unsigned y = 0; y < height; y++) {
293 for(unsigned x = 0; x < width; x++) {
294 dst_pixel[x] *= src_pixel[x];
295 }
296 // next line
297 src_pixel += pitch;
298 dst_pixel += pitch;
299 }
300
301 FreeImage_Unload(L);
302
303 // PHI(k+1) is no longer needed
304 FreeImage_Unload(phi[k+1]);
305 phi[k+1] = NULL;
306 }
307
308 // next level
309 }
310
311 // get the final result and return
312 FIBITMAP *dst = phi[0];
313
314 free(phi);
315
316 return dst;
317
318 } catch(int) {
319 if(phi) {
320 for(int k = nlevels-1; k >= 0; k--) {
321 if(phi[k]) FreeImage_Unload(phi[k]);
322 }
323 free(phi);
324 }
325 return NULL;
326 }
327 }
328
329 /**
330 Compute gradients in x and y directions, attenuate them with the attenuation matrix,
331 then compute the divergence div G from the attenuated gradient.
332 @param H Normalized luminance
333 @param PHI Attenuation matrix
334 @return Returns the divergence matrix if successful, returns NULL otherwise
335 */
Divergence(FIBITMAP * H,FIBITMAP * PHI)336 static FIBITMAP* Divergence(FIBITMAP *H, FIBITMAP *PHI) {
337 FIBITMAP *Gx = NULL, *Gy = NULL, *divG = NULL;
338 float *phi, *h, *gx, *gy, *divg;
339
340 try {
341 const FREE_IMAGE_TYPE image_type = FreeImage_GetImageType(H);
342 if(image_type != FIT_FLOAT) throw(1);
343
344 const unsigned width = FreeImage_GetWidth(H);
345 const unsigned height = FreeImage_GetHeight(H);
346
347 Gx = FreeImage_AllocateT(image_type, width, height);
348 if(!Gx) throw(1);
349 Gy = FreeImage_AllocateT(image_type, width, height);
350 if(!Gy) throw(1);
351
352 const unsigned pitch = FreeImage_GetPitch(H) / sizeof(float);
353
354 // perform gradient attenuation
355
356 phi = (float*)FreeImage_GetBits(PHI);
357 h = (float*)FreeImage_GetBits(H);
358 gx = (float*)FreeImage_GetBits(Gx);
359 gy = (float*)FreeImage_GetBits(Gy);
360
361 for(unsigned y = 0; y < height; y++) {
362 const unsigned s = (y+1 == height ? y : y+1);
363 for(unsigned x = 0; x < width; x++) {
364 const unsigned e = (x+1 == width ? x : x+1);
365 // forward difference
366 const unsigned index = y*pitch + x;
367 const float phi_xy = phi[index];
368 const float h_xy = h[index];
369 gx[x] = (h[y*pitch+e] - h_xy) * phi_xy; // [H(x+1, y) - H(x, y)] * PHI(x, y)
370 gy[x] = (h[s*pitch+x] - h_xy) * phi_xy; // [H(x, y+1) - H(x, y)] * PHI(x, y)
371 }
372 // next line
373 gx += pitch;
374 gy += pitch;
375 }
376
377 // calculate the divergence
378
379 divG = FreeImage_AllocateT(image_type, width, height);
380 if(!divG) throw(1);
381
382 gx = (float*)FreeImage_GetBits(Gx);
383 gy = (float*)FreeImage_GetBits(Gy);
384 divg = (float*)FreeImage_GetBits(divG);
385
386 for(unsigned y = 0; y < height; y++) {
387 for(unsigned x = 0; x < width; x++) {
388 // backward difference approximation
389 // divG = Gx(x, y) - Gx(x-1, y) + Gy(x, y) - Gy(x, y-1)
390 const unsigned index = y*pitch + x;
391 divg[index] = gx[index] + gy[index];
392 if(x > 0) divg[index] -= gx[index-1];
393 if(y > 0) divg[index] -= gy[index-pitch];
394 }
395 }
396
397 // no longer needed ...
398 FreeImage_Unload(Gx);
399 FreeImage_Unload(Gy);
400
401 // return the divergence
402 return divG;
403
404 } catch(int) {
405 if(Gx) FreeImage_Unload(Gx);
406 if(Gy) FreeImage_Unload(Gy);
407 if(divG) FreeImage_Unload(divG);
408 return NULL;
409 }
410 }
411
412 /**
413 Given the luminance channel, find max & min luminance values,
414 normalize to range 0..100 and take the logarithm.
415 @param Y Image luminance
416 @return Returns the normalized luminance H if successful, returns NULL otherwise
417 */
LogLuminance(FIBITMAP * Y)418 static FIBITMAP* LogLuminance(FIBITMAP *Y) {
419 FIBITMAP *H = NULL;
420
421 try {
422 // get the luminance channel
423 FIBITMAP *H = FreeImage_Clone(Y);
424 if(!H) throw(1);
425
426 const unsigned width = FreeImage_GetWidth(H);
427 const unsigned height = FreeImage_GetHeight(H);
428 const unsigned pitch = FreeImage_GetPitch(H);
429
430 // find max & min luminance values
431 float maxLum = -1e20F, minLum = 1e20F;
432
433 BYTE *bits = (BYTE*)FreeImage_GetBits(H);
434 for(unsigned y = 0; y < height; y++) {
435 const float *pixel = (float*)bits;
436 for(unsigned x = 0; x < width; x++) {
437 const float value = pixel[x];
438 maxLum = (maxLum < value) ? value : maxLum; // max Luminance in the scene
439 minLum = (minLum < value) ? minLum : value; // min Luminance in the scene
440 }
441 // next line
442 bits += pitch;
443 }
444 if(maxLum == minLum) throw(1);
445
446 // normalize to range 0..100 and take the logarithm
447 const float scale = 100.F / (maxLum - minLum);
448 bits = (BYTE*)FreeImage_GetBits(H);
449 for(unsigned y = 0; y < height; y++) {
450 float *pixel = (float*)bits;
451 for(unsigned x = 0; x < width; x++) {
452 const float value = (pixel[x] - minLum) * scale;
453 pixel[x] = log(value + EPSILON);
454 }
455 // next line
456 bits += pitch;
457 }
458
459 return H;
460
461 } catch(int) {
462 if(H) FreeImage_Unload(H);
463 return NULL;
464 }
465 }
466
467 /**
468 Given a normalized luminance, perform exponentiation and recover the log compressed image
469 @param Y Input/Output luminance image
470 */
ExpLuminance(FIBITMAP * Y)471 static void ExpLuminance(FIBITMAP *Y) {
472 const unsigned width = FreeImage_GetWidth(Y);
473 const unsigned height = FreeImage_GetHeight(Y);
474 const unsigned pitch = FreeImage_GetPitch(Y);
475
476 BYTE *bits = (BYTE*)FreeImage_GetBits(Y);
477 for(unsigned y = 0; y < height; y++) {
478 float *pixel = (float*)bits;
479 for(unsigned x = 0; x < width; x++) {
480 pixel[x] = exp(pixel[x]) - EPSILON;
481 }
482 bits += pitch;
483 }
484 }
485
486 // --------------------------------------------------------------------------
487
488 /**
489 Gradient Domain HDR tone mapping operator
490 @param Y Image luminance values
491 @param alpha Parameter alpha of the paper (suggested value is 0.1)
492 @param beta Parameter beta of the paper (suggested value is between 0.8 and 0.9)
493 @return returns the tone mapped luminance
494 */
tmoFattal02(FIBITMAP * Y,float alpha,float beta)495 static FIBITMAP* tmoFattal02(FIBITMAP *Y, float alpha, float beta) {
496 const unsigned MIN_PYRAMID_SIZE = 32; // minimun size (width or height) of the coarsest level of the pyramid
497
498 FIBITMAP *H = NULL;
499 FIBITMAP **pyramid = NULL;
500 FIBITMAP **gradients = NULL;
501 FIBITMAP *phy = NULL;
502 FIBITMAP *divG = NULL;
503 FIBITMAP *U = NULL;
504 float *avgGrad = NULL;
505
506 int k;
507 int nlevels = 0;
508
509 try {
510 // get the normalized luminance
511 FIBITMAP *H = LogLuminance(Y);
512 if(!H) throw(1);
513
514 // get the number of levels for the pyramid
515 const unsigned width = FreeImage_GetWidth(H);
516 const unsigned height = FreeImage_GetHeight(H);
517 unsigned minsize = MIN(width, height);
518 while(minsize >= MIN_PYRAMID_SIZE) {
519 nlevels++;
520 minsize /= 2;
521 }
522
523 // create the Gaussian pyramid
524 pyramid = (FIBITMAP**)malloc(nlevels * sizeof(FIBITMAP*));
525 if(!pyramid) throw(1);
526 memset(pyramid, 0, nlevels * sizeof(FIBITMAP*));
527
528 if(!GaussianPyramid(H, pyramid, nlevels)) throw(1);
529
530 // calculate gradient magnitude and its average value on each pyramid level
531 gradients = (FIBITMAP**)malloc(nlevels * sizeof(FIBITMAP*));
532 if(!gradients) throw(1);
533 memset(gradients, 0, nlevels * sizeof(FIBITMAP*));
534 avgGrad = (float*)malloc(nlevels * sizeof(float));
535 if(!avgGrad) throw(1);
536
537 if(!GradientPyramid(pyramid, nlevels, gradients, avgGrad)) throw(1);
538
539 // free the Gaussian pyramid
540 for(k = 0; k < nlevels; k++) {
541 if(pyramid[k]) FreeImage_Unload(pyramid[k]);
542 }
543 free(pyramid); pyramid = NULL;
544
545 // compute the gradient attenuation function PHI(x, y)
546 phy = PhiMatrix(gradients, avgGrad, nlevels, alpha, beta);
547 if(!phy) throw(1);
548
549 // free the gradient pyramid
550 for(k = 0; k < nlevels; k++) {
551 if(gradients[k]) FreeImage_Unload(gradients[k]);
552 }
553 free(gradients); gradients = NULL;
554 free(avgGrad); avgGrad = NULL;
555
556 // compute gradients in x and y directions, attenuate them with the attenuation matrix,
557 // then compute the divergence div G from the attenuated gradient.
558 divG = Divergence(H, phy);
559 if(!divG) throw(1);
560
561 // H & phy no longer needed
562 FreeImage_Unload(H); H = NULL;
563 FreeImage_Unload(phy); phy = NULL;
564
565 // solve the PDE (Poisson equation) using a multigrid solver and 3 cycles
566 FIBITMAP *U = FreeImage_MultigridPoissonSolver(divG, 3);
567 if(!U) throw(1);
568
569 FreeImage_Unload(divG);
570
571 // perform exponentiation and recover the log compressed image
572 ExpLuminance(U);
573
574 return U;
575
576 } catch(int) {
577 if(H) FreeImage_Unload(H);
578 if(pyramid) {
579 for(int i = 0; i < nlevels; i++) {
580 if(pyramid[i]) FreeImage_Unload(pyramid[i]);
581 }
582 free(pyramid);
583 }
584 if(gradients) {
585 for(int i = 0; i < nlevels; i++) {
586 if(gradients[i]) FreeImage_Unload(gradients[i]);
587 }
588 free(gradients);
589 }
590 if(avgGrad) free(avgGrad);
591 if(phy) FreeImage_Unload(phy);
592 if(divG) FreeImage_Unload(divG);
593 if(U) FreeImage_Unload(U);
594
595 return NULL;
596 }
597 }
598
599 // ----------------------------------------------------------
600 // Main algorithm
601 // ----------------------------------------------------------
602
603 /**
604 Apply the Gradient Domain High Dynamic Range Compression to a RGBF image and convert to 24-bit RGB
605 @param dib Input RGBF / RGB16 image
606 @param color_saturation Color saturation (s parameter in the paper) in [0.4..0.6]
607 @param attenuation Atenuation factor (beta parameter in the paper) in [0.8..0.9]
608 @return Returns a 24-bit RGB image if successful, returns NULL otherwise
609 */
610 FIBITMAP* DLL_CALLCONV
FreeImage_TmoFattal02(FIBITMAP * dib,double color_saturation,double attenuation)611 FreeImage_TmoFattal02(FIBITMAP *dib, double color_saturation, double attenuation) {
612 const float alpha = 0.1F; // parameter alpha = 0.1
613 const float beta = (float)MAX(0.8, MIN(0.9, attenuation)); // parameter beta = [0.8..0.9]
614 const float s = (float)MAX(0.4, MIN(0.6, color_saturation));// exponent s controls color saturation = [0.4..0.6]
615
616 FIBITMAP *src = NULL;
617 FIBITMAP *Yin = NULL;
618 FIBITMAP *Yout = NULL;
619 FIBITMAP *dst = NULL;
620
621 if(!FreeImage_HasPixels(dib)) return NULL;
622
623 try {
624
625 // convert to RGBF
626 src = FreeImage_ConvertToRGBF(dib);
627 if(!src) throw(1);
628
629 // get the luminance channel
630 Yin = ConvertRGBFToY(src);
631 if(!Yin) throw(1);
632
633 // perform the tone mapping
634 Yout = tmoFattal02(Yin, alpha, beta);
635 if(!Yout) throw(1);
636
637 // clip low and high values and normalize to [0..1]
638 //NormalizeY(Yout, 0.001F, 0.995F);
639 NormalizeY(Yout, 0, 1);
640
641 // compress the dynamic range
642
643 const unsigned width = FreeImage_GetWidth(src);
644 const unsigned height = FreeImage_GetHeight(src);
645
646 const unsigned rgb_pitch = FreeImage_GetPitch(src);
647 const unsigned y_pitch = FreeImage_GetPitch(Yin);
648
649 BYTE *bits = (BYTE*)FreeImage_GetBits(src);
650 BYTE *bits_yin = (BYTE*)FreeImage_GetBits(Yin);
651 BYTE *bits_yout = (BYTE*)FreeImage_GetBits(Yout);
652
653 for(unsigned y = 0; y < height; y++) {
654 float *Lin = (float*)bits_yin;
655 float *Lout = (float*)bits_yout;
656 float *color = (float*)bits;
657 for(unsigned x = 0; x < width; x++) {
658 for(unsigned c = 0; c < 3; c++) {
659 *color = (Lin[x] > 0) ? pow(*color/Lin[x], s) * Lout[x] : 0;
660 color++;
661 }
662 }
663 bits += rgb_pitch;
664 bits_yin += y_pitch;
665 bits_yout += y_pitch;
666 }
667
668 // not needed anymore
669 FreeImage_Unload(Yin); Yin = NULL;
670 FreeImage_Unload(Yout); Yout = NULL;
671
672 // clamp image highest values to display white, then convert to 24-bit RGB
673 dst = ClampConvertRGBFTo24(src);
674
675 // clean-up and return
676 FreeImage_Unload(src); src = NULL;
677
678 // copy metadata from src to dst
679 FreeImage_CloneMetadata(dst, dib);
680
681 return dst;
682
683 } catch(int) {
684 if(src) FreeImage_Unload(src);
685 if(Yin) FreeImage_Unload(Yin);
686 if(Yout) FreeImage_Unload(Yout);
687 return NULL;
688 }
689 }
690