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 /**
40 nominal CRT primaries
41 */
42 /*
43 static const float CIE_x_r = 0.640F;
44 static const float CIE_y_r = 0.330F;
45 static const float CIE_x_g = 0.290F;
46 static const float CIE_y_g = 0.600F;
47 static const float CIE_x_b = 0.150F;
48 static const float CIE_y_b = 0.060F;
49 static const float CIE_x_w = 0.3333F; // use true white
50 static const float CIE_y_w = 0.3333F;
51 */
52 /**
53 sRGB primaries
54 */
55 static const float CIE_x_r = 0.640F;
56 static const float CIE_y_r = 0.330F;
57 static const float CIE_x_g = 0.300F;
58 static const float CIE_y_g = 0.600F;
59 static const float CIE_x_b = 0.150F;
60 static const float CIE_y_b = 0.060F;
61 static const float CIE_x_w = 0.3127F; // Illuminant D65
62 static const float CIE_y_w = 0.3290F;
63
64 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) );
65 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) );
66 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) );
67 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) );
68
69 /**
70 RGB to XYZ (no white balance)
71 */
72 static 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 /**
88 XYZ to RGB (no white balance)
89 */
90 static 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 /**
106 This gives approximately the following matrices :
107
108 static 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 };
113 static 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
122 static const float EPSILON = 1e-06F;
123 static const float INF = 1e+10F;
124
125 /**
126 Convert in-place floating point RGB data to Yxy.<br>
127 On 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 */
131 BOOL
ConvertInPlaceRGBFToYxy(FIBITMAP * dib)132 ConvertInPlaceRGBFToYxy(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 /**
170 Convert in-place Yxy image to floating point RGB data.<br>
171 On 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 */
175 BOOL
ConvertInPlaceYxyToRGBF(FIBITMAP * dib)176 ConvertInPlaceYxyToRGBF(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 /**
221 Get the maximum, minimum and average luminance.<br>
222 On 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 */
229 BOOL
LuminanceFromYxy(FIBITMAP * Yxy,float * maxLum,float * minLum,float * worldLum)230 LuminanceFromYxy(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 = MAX(0.0F, pixel[x].red);// avoid negative values
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 /**
266 Clamp RGBF image highest values to display white,
267 then convert to 24-bit RGB
268 */
269 FIBITMAP*
ClampConvertRGBFTo24(FIBITMAP * src)270 ClampConvertRGBFTo24(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 /**
307 Extract the luminance channel L from a RGBF image.
308 Luminance is calculated from the sRGB model (RGB2XYZ matrix)
309 using a D65 white point :
310 L = ( 0.2126 * r ) + ( 0.7152 * g ) + ( 0.0722 * b )
311 Reference :
312 A Standard Default Color Space for the Internet - sRGB.
313 [online] http://www.w3.org/Graphics/Color/sRGB
314 */
315 FIBITMAP*
ConvertRGBFToY(FIBITMAP * src)316 ConvertRGBFToY(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 /**
349 Get 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 */
358 BOOL
LuminanceFromY(FIBITMAP * dib,float * maxLum,float * minLum,float * Lav,float * Llav)359 LuminanceFromY(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
findMaxMinPercentile(FIBITMAP * Y,float minPrct,float * minLum,float maxPrct,float * maxLum)397 static 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 /**
423 Clipping function<br>
424 Remove any extremely bright and/or extremely dark pixels
425 and normalize between 0 and 1.
426 @param Y Input/Output image
427 @param minPrct Minimum percentile
428 @param maxPrct Maximum percentile
429 */
430 void
NormalizeY(FIBITMAP * Y,float minPrct,float maxPrct)431 NormalizeY(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 }
480