PageRenderTime 39ms CodeModel.GetById 1ms app.highlight 33ms RepoModel.GetById 1ms app.codeStats 1ms

/src/FreeImage/Source/FreeImage/tmoColorConvert.cpp

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