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
22 #include "FreeImage.h"
23 #include "Utilities.h"
24 #include "ToneMapping.h"
25
26 // ----------------------------------------------------------
27 // Logarithmic mapping operator
28 // Reference:
29 // [1] F. Drago, K. Myszkowski, T. Annen, and N. Chiba,
30 // Adaptive Logarithmic Mapping for Displaying High Contrast Scenes,
31 // Eurographics 2003.
32 // ----------------------------------------------------------
33
34 /**
35 Bias function
36 */
37 static inline double
biasFunction(const double b,const double x)38 biasFunction(const double b, const double x) {
39 return pow (x, b); // pow(x, log(bias)/log(0.5)
40 }
41
42 /**
43 Pad� approximation of log(x + 1)
44 x(6+x)/(6+4x) good if x < 1
45 x*(6 + 0.7662x)/(5.9897 + 3.7658x) between 1 and 2
46 See http://www.nezumi.demon.co.uk/consult/logx.htm
47 */
48 static inline double
pade_log(const double x)49 pade_log(const double x) {
50 if(x < 1) {
51 return (x * (6 + x) / (6 + 4 * x));
52 } else if(x < 2) {
53 return (x * (6 + 0.7662 * x) / (5.9897 + 3.7658 * x));
54 }
55 return log(x + 1);
56 }
57
58 /**
59 Log mapping operator
60 @param dib Input / Output Yxy image
61 @param maxLum Maximum luminance
62 @param avgLum Average luminance (world adaptation luminance)
63 @param biasParam Bias parameter (a zero value default to 0.85)
64 @param exposure Exposure parameter (default to 0)
65 @return Returns TRUE if successful, returns FALSE otherwise
66 */
67 static BOOL
ToneMappingDrago03(FIBITMAP * dib,const float maxLum,const float avgLum,float biasParam,const float exposure)68 ToneMappingDrago03(FIBITMAP *dib, const float maxLum, const float avgLum, float biasParam, const float exposure) {
69 const float LOG05 = -0.693147F; // log(0.5)
70
71 double Lmax, divider, interpol, biasP;
72 unsigned x, y;
73 double L;
74
75 if(FreeImage_GetImageType(dib) != FIT_RGBF)
76 return FALSE;
77
78 const unsigned width = FreeImage_GetWidth(dib);
79 const unsigned height = FreeImage_GetHeight(dib);
80 const unsigned pitch = FreeImage_GetPitch(dib);
81
82
83 // arbitrary Bias Parameter
84 if(biasParam == 0)
85 biasParam = 0.85F;
86
87 // normalize maximum luminance by average luminance
88 Lmax = maxLum / avgLum;
89
90 divider = log10(Lmax+1);
91 biasP = log(biasParam)/LOG05;
92
93 #if !defined(DRAGO03_FAST)
94
95 /**
96 Normal tone mapping of every pixel
97 further acceleration is obtained by a Pad� approximation of log(x + 1)
98 */
99 BYTE *bits = (BYTE*)FreeImage_GetBits(dib);
100 for(y = 0; y < height; y++) {
101 FIRGBF *pixel = (FIRGBF*)bits;
102 for(x = 0; x < width; x++) {
103 double Yw = pixel[x].red / avgLum;
104 Yw *= exposure;
105 interpol = log(2 + biasFunction(biasP, Yw / Lmax) * 8);
106 L = pade_log(Yw);// log(Yw + 1)
107 pixel[x].red = (float)((L / interpol) / divider);
108 }
109 // next line
110 bits += pitch;
111 }
112
113 #else
114 unsigned index;
115 int i, j;
116
117 unsigned max_width = width - (width % 3);
118 unsigned max_height = height - (height % 3);
119 unsigned fpitch = pitch / sizeof(FIRGBF);
120
121 /**
122 fast tone mapping
123 split the image into 3x3 pixel tiles and perform the computation for each group of 9 pixels
124 further acceleration is obtained by a Pad� approximation of log(x + 1)
125 => produce artifacts and not so faster, so the code has been disabled
126 */
127 #define PIXEL(x, y) image[y*fpitch + x].red
128
129 FIRGBF *image = (FIRGBF*)FreeImage_GetBits(dib);
130 for(y = 0; y < max_height; y += 3) {
131 for(x = 0; x < max_width; x += 3) {
132 double average = 0;
133 for(i = 0; i < 3; i++) {
134 for(j = 0; j < 3; j++) {
135 index = (y + i)*fpitch + (x + j);
136 image[index].red /= (float)avgLum;
137 image[index].red *= exposure;
138 average += image[index].red;
139 }
140 }
141 average = average / 9 - PIXEL(x, y);
142 if(average > -1 && average < 1) {
143 interpol = log(2 + pow(PIXEL(x + 1, y + 1) / Lmax, biasP) * 8);
144 for(i = 0; i < 3; i++) {
145 for(j = 0; j < 3; j++) {
146 index = (y + i)*fpitch + (x + j);
147 L = pade_log(image[index].red);// log(image[index].red + 1)
148 image[index].red = (float)((L / interpol) / divider);
149 }
150 }
151 }
152 else {
153 for(i = 0; i < 3; i++) {
154 for(j = 0; j < 3; j++) {
155 index = (y + i)*fpitch + (x + j);
156 interpol = log(2 + pow(image[index].red / Lmax, biasP) * 8);
157 L = pade_log(image[index].red);// log(image[index].red + 1)
158 image[index].red = (float)((L / interpol) / divider);
159 }
160 }
161 }
162 } //x
163 } // y
164
165 /**
166 Normal tone mapping of every pixel for the remaining right and bottom bands
167 */
168 BYTE *bits;
169
170 // right band
171 bits = (BYTE*)FreeImage_GetBits(dib);
172 for(y = 0; y < height; y++) {
173 FIRGBF *pixel = (FIRGBF*)bits;
174 for(x = max_width; x < width; x++) {
175 double Yw = pixel[x].red / avgLum;
176 Yw *= exposure;
177 interpol = log(2 + biasFunction(biasP, Yw / Lmax) * 8);
178 L = pade_log(Yw);// log(Yw + 1)
179 pixel[x].red = (float)((L / interpol) / divider);
180 }
181 // next line
182 bits += pitch;
183 }
184 // bottom band
185 bits = (BYTE*)FreeImage_GetBits(dib);
186 for(y = max_height; y < height; y++) {
187 FIRGBF *pixel = (FIRGBF*)bits;
188 for(x = 0; x < max_width; x++) {
189 double Yw = pixel[x].red / avgLum;
190 Yw *= exposure;
191 interpol = log(2 + biasFunction(biasP, Yw / Lmax) * 8);
192 L = pade_log(Yw);// log(Yw + 1)
193 pixel[x].red = (float)((L / interpol) / divider);
194 }
195 // next line
196 bits += pitch;
197 }
198
199 #endif // DRAGO03_FAST
200
201 return TRUE;
202 }
203
204 /**
205 Custom gamma correction based on the ITU-R BT.709 standard
206 @param dib RGBF image to be corrected
207 @param gammaval Gamma value (2.2 is a good default value)
208 @return Returns TRUE if successful, returns FALSE otherwise
209 */
210 static BOOL
REC709GammaCorrection(FIBITMAP * dib,const float gammaval)211 REC709GammaCorrection(FIBITMAP *dib, const float gammaval) {
212 if(FreeImage_GetImageType(dib) != FIT_RGBF)
213 return FALSE;
214
215 float slope = 4.5F;
216 float start = 0.018F;
217
218 const float fgamma = (float)((0.45 / gammaval) * 2);
219 if(gammaval >= 2.1F) {
220 start = (float)(0.018 / ((gammaval - 2) * 7.5));
221 slope = (float)(4.5 * ((gammaval - 2) * 7.5));
222 } else if (gammaval <= 1.9F) {
223 start = (float)(0.018 * ((2 - gammaval) * 7.5));
224 slope = (float)(4.5 / ((2 - gammaval) * 7.5));
225 }
226
227 const unsigned width = FreeImage_GetWidth(dib);
228 const unsigned height = FreeImage_GetHeight(dib);
229 const unsigned pitch = FreeImage_GetPitch(dib);
230
231 BYTE *bits = (BYTE*)FreeImage_GetBits(dib);
232 for(unsigned y = 0; y < height; y++) {
233 float *pixel = (float*)bits;
234 for(unsigned x = 0; x < width; x++) {
235 for(int i = 0; i < 3; i++) {
236 *pixel = (*pixel <= start) ? *pixel * slope : (1.099F * pow(*pixel, fgamma) - 0.099F);
237 pixel++;
238 }
239 }
240 bits += pitch;
241 }
242
243 return TRUE;
244 }
245
246 // ----------------------------------------------------------
247 // Main algorithm
248 // ----------------------------------------------------------
249
250 /**
251 Apply the Adaptive Logarithmic Mapping operator to a HDR image and convert to 24-bit RGB
252 @param src Input RGB16 or RGB[A]F image
253 @param gamma Gamma correction (gamma > 0). 1 means no correction, 2.2 in the original paper.
254 @param exposure Exposure parameter (0 means no correction, 0 in the original paper)
255 @return Returns a 24-bit RGB image if successful, returns NULL otherwise
256 */
257 FIBITMAP* DLL_CALLCONV
FreeImage_TmoDrago03(FIBITMAP * src,double gamma,double exposure)258 FreeImage_TmoDrago03(FIBITMAP *src, double gamma, double exposure) {
259 float maxLum, minLum, avgLum;
260
261 if(!FreeImage_HasPixels(src)) return NULL;
262
263 // working RGBF variable
264 FIBITMAP *dib = NULL;
265
266 dib = FreeImage_ConvertToRGBF(src);
267 if(!dib) return NULL;
268
269 // default algorithm parameters
270 const float biasParam = 0.85F;
271 const float expoParam = (float)pow(2.0, exposure); //default exposure is 1, 2^0
272
273 // convert to Yxy
274 ConvertInPlaceRGBFToYxy(dib);
275 // get the luminance
276 LuminanceFromYxy(dib, &maxLum, &minLum, &avgLum);
277 // perform the tone mapping
278 ToneMappingDrago03(dib, maxLum, avgLum, biasParam, expoParam);
279 // convert back to RGBF
280 ConvertInPlaceYxyToRGBF(dib);
281 if(gamma != 1) {
282 // perform gamma correction
283 REC709GammaCorrection(dib, (float)gamma);
284 }
285 // clamp image highest values to display white, then convert to 24-bit RGB
286 FIBITMAP *dst = ClampConvertRGBFTo24(dib);
287
288 // clean-up and return
289 FreeImage_Unload(dib);
290
291 // copy metadata from src to dst
292 FreeImage_CloneMetadata(dst, src);
293
294 return dst;
295 }
296