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

  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/**
 35Bias function
 36*/
 37static inline double 
 38biasFunction(const double b, const double x) {
 39	return pow (x, b);		// pow(x, log(bias)/log(0.5)
 40}
 41
 42/**
 43Padé approximation of log(x + 1)
 44x(6+x)/(6+4x) good if x < 1
 45x*(6 + 0.7662x)/(5.9897 + 3.7658x) between 1 and 2
 46See http://www.nezumi.demon.co.uk/consult/logx.htm
 47*/
 48static inline double 
 49pade_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/**
 59Log 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*/
 67static BOOL 
 68ToneMappingDrago03(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/**
205Custom 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*/
210static BOOL 
211REC709GammaCorrection(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/**
251Apply 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*/
257FIBITMAP* DLL_CALLCONV 
258FreeImage_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}