From 87d155ba9f4e699293adc65954c89f20fc65b36d Mon Sep 17 00:00:00 2001 From: Kirill Volinsky Date: Sat, 22 Dec 2012 16:38:58 +0000 Subject: AdvaImg plugin folder renaming git-svn-id: http://svn.miranda-ng.org/main/trunk@2801 1316c22d-e87f-b044-9b9b-93d7a3e3ba9c --- plugins/AdvaImg/src/FreeImage/tmoFattal02.cpp | 689 ++++++++++++++++++++++++++ 1 file changed, 689 insertions(+) create mode 100644 plugins/AdvaImg/src/FreeImage/tmoFattal02.cpp (limited to 'plugins/AdvaImg/src/FreeImage/tmoFattal02.cpp') diff --git a/plugins/AdvaImg/src/FreeImage/tmoFattal02.cpp b/plugins/AdvaImg/src/FreeImage/tmoFattal02.cpp new file mode 100644 index 0000000000..bcd4524523 --- /dev/null +++ b/plugins/AdvaImg/src/FreeImage/tmoFattal02.cpp @@ -0,0 +1,689 @@ +// ========================================================== +// Tone mapping operator (Fattal, 2002) +// +// Design and implementation by +// - Hervé Drolon (drolon@infonie.fr) +// +// This file is part of FreeImage 3 +// +// COVERED CODE IS PROVIDED UNDER THIS LICENSE ON AN "AS IS" BASIS, WITHOUT WARRANTY +// OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, WITHOUT LIMITATION, WARRANTIES +// THAT THE COVERED CODE IS FREE OF DEFECTS, MERCHANTABLE, FIT FOR A PARTICULAR PURPOSE +// OR NON-INFRINGING. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE COVERED +// CODE IS WITH YOU. SHOULD ANY COVERED CODE PROVE DEFECTIVE IN ANY RESPECT, YOU (NOT +// THE INITIAL DEVELOPER OR ANY OTHER CONTRIBUTOR) ASSUME THE COST OF ANY NECESSARY +// SERVICING, REPAIR OR CORRECTION. THIS DISCLAIMER OF WARRANTY CONSTITUTES AN ESSENTIAL +// PART OF THIS LICENSE. NO USE OF ANY COVERED CODE IS AUTHORIZED HEREUNDER EXCEPT UNDER +// THIS DISCLAIMER. +// +// Use at your own risk! +// ========================================================== + +#include "FreeImage.h" +#include "Utilities.h" +#include "ToneMapping.h" + +// ---------------------------------------------------------- +// Gradient domain HDR compression +// Reference: +// [1] R. Fattal, D. Lischinski, and M.Werman, +// Gradient domain high dynamic range compression, +// ACM Transactions on Graphics, special issue on Proc. of ACM SIGGRAPH 2002, +// San Antonio, Texas, vol. 21(3), pp. 257-266, 2002. +// ---------------------------------------------------------- + +static const float EPSILON = 1e-4F; + +/** +Performs a 5 by 5 gaussian filtering using two 1D convolutions, +followed by a subsampling by 2. +@param dib Input image +@return Returns a blurred image of size SIZE(dib)/2 +@see GaussianPyramid +*/ +static FIBITMAP* GaussianLevel5x5(FIBITMAP *dib) { + FIBITMAP *h_dib = NULL, *v_dib = NULL, *dst = NULL; + float *src_pixel, *dst_pixel; + + try { + const FREE_IMAGE_TYPE image_type = FreeImage_GetImageType(dib); + if(image_type != FIT_FLOAT) throw(1); + + const unsigned width = FreeImage_GetWidth(dib); + const unsigned height = FreeImage_GetHeight(dib); + + h_dib = FreeImage_AllocateT(image_type, width, height); + v_dib = FreeImage_AllocateT(image_type, width, height); + if (!h_dib || !v_dib) throw(1); + + const unsigned pitch = FreeImage_GetPitch(dib) / sizeof(float); + + // horizontal convolution dib -> h_dib + + src_pixel = (float*)FreeImage_GetBits(dib); + dst_pixel = (float*)FreeImage_GetBits(h_dib); + + for(unsigned y = 0; y < height; y++) { + // work on line y + for(unsigned x = 2; x < width - 2; x++) { + dst_pixel[x] = src_pixel[x-2] + src_pixel[x+2] + 4 * (src_pixel[x-1] + src_pixel[x+1]) + 6 * src_pixel[x]; + dst_pixel[x] /= 16; + } + // boundary mirroring + dst_pixel[0] = (2 * src_pixel[2] + 8 * src_pixel[1] + 6 * src_pixel[0]) / 16; + dst_pixel[1] = (src_pixel[3] + 4 * (src_pixel[0] + src_pixel[2]) + 7 * src_pixel[1]) / 16; + dst_pixel[width-2] = (src_pixel[width-4] + 5 * src_pixel[width-1] + 4 * src_pixel[width-3] + 6 * src_pixel[width-2]) / 16; + dst_pixel[width-1] = (src_pixel[width-3] + 5 * src_pixel[width-2] + 10 * src_pixel[width-1]) / 16; + + // next line + src_pixel += pitch; + dst_pixel += pitch; + } + + // vertical convolution h_dib -> v_dib + + src_pixel = (float*)FreeImage_GetBits(h_dib); + dst_pixel = (float*)FreeImage_GetBits(v_dib); + + for(unsigned x = 0; x < width; x++) { + // work on column x + for(unsigned y = 2; y < height - 2; y++) { + const unsigned index = y*pitch + x; + 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]; + dst_pixel[index] /= 16; + } + // boundary mirroring + dst_pixel[x] = (2 * src_pixel[x+2*pitch] + 8 * src_pixel[x+pitch] + 6 * src_pixel[x]) / 16; + dst_pixel[x+pitch] = (src_pixel[x+3*pitch] + 4 * (src_pixel[x] + src_pixel[x+2*pitch]) + 7 * src_pixel[x+pitch]) / 16; + 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; + 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; + } + + FreeImage_Unload(h_dib); h_dib = NULL; + + // perform downsampling + + dst = FreeImage_Rescale(v_dib, width/2, height/2, FILTER_BILINEAR); + + FreeImage_Unload(v_dib); + + return dst; + + } catch(int) { + if(h_dib) FreeImage_Unload(h_dib); + if(v_dib) FreeImage_Unload(v_dib); + if(dst) FreeImage_Unload(dst); + return NULL; + } +} + +/** +Compute a Gaussian pyramid using the specified number of levels. +@param H Original bitmap +@param pyramid Resulting pyramid array +@param nlevels Number of resolution levels +@return Returns TRUE if successful, returns FALSE otherwise +*/ +static BOOL GaussianPyramid(FIBITMAP *H, FIBITMAP **pyramid, int nlevels) { + try { + // first level is the original image + pyramid[0] = FreeImage_Clone(H); + if(pyramid[0] == NULL) throw(1); + // compute next levels + for(int k = 1; k < nlevels; k++) { + pyramid[k] = GaussianLevel5x5(pyramid[k-1]); + if(pyramid[k] == NULL) throw(1); + } + return TRUE; + } catch(int) { + for(int k = 0; k < nlevels; k++) { + if(pyramid[k] != NULL) { + FreeImage_Unload(pyramid[k]); + pyramid[k] = NULL; + } + } + return FALSE; + } +} + +/** +Compute the gradient magnitude of an input image H using central differences, +and returns the average gradient. +@param H Input image +@param avgGrad [out] Average gradient +@param k Level number +@return Returns the gradient magnitude if successful, returns NULL otherwise +@see GradientPyramid +*/ +static FIBITMAP* GradientLevel(FIBITMAP *H, float *avgGrad, int k) { + FIBITMAP *G = NULL; + + try { + const FREE_IMAGE_TYPE image_type = FreeImage_GetImageType(H); + if(image_type != FIT_FLOAT) throw(1); + + const unsigned width = FreeImage_GetWidth(H); + const unsigned height = FreeImage_GetHeight(H); + + G = FreeImage_AllocateT(image_type, width, height); + if (!G) throw(1); + + const unsigned pitch = FreeImage_GetPitch(H) / sizeof(float); + + const float divider = (float)(1 << (k + 1)); + float average = 0; + + float *src_pixel = (float*)FreeImage_GetBits(H); + float *dst_pixel = (float*)FreeImage_GetBits(G); + + for(unsigned y = 0; y < height; y++) { + const unsigned n = (y == 0 ? 0 : y-1); + const unsigned s = (y+1 == height ? y : y+1); + for(unsigned x = 0; x < width; x++) { + const unsigned w = (x == 0 ? 0 : x-1); + const unsigned e = (x+1 == width ? x : x+1); + // central difference + 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) + 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) + // gradient + dst_pixel[x] = sqrt(gx*gx + gy*gy); + // average gradient + average += dst_pixel[x]; + } + // next line + dst_pixel += pitch; + } + + *avgGrad = average / (width * height); + + return G; + + } catch(int) { + if(G) FreeImage_Unload(G); + return NULL; + } +} + +/** +Calculate gradient magnitude and its average value on each pyramid level +@param pyramid Gaussian pyramid (nlevels levels) +@param nlevels Number of levels +@param gradients [out] Gradient pyramid (nlevels levels) +@param avgGrad [out] Average gradient on each level (array of size nlevels) +@return Returns TRUE if successful, returns FALSE otherwise +*/ +static BOOL GradientPyramid(FIBITMAP **pyramid, int nlevels, FIBITMAP **gradients, float *avgGrad) { + try { + for(int k = 0; k < nlevels; k++) { + FIBITMAP *Hk = pyramid[k]; + gradients[k] = GradientLevel(Hk, &avgGrad[k], k); + if(gradients[k] == NULL) throw(1); + } + return TRUE; + } catch(int) { + for(int k = 0; k < nlevels; k++) { + if(gradients[k] != NULL) { + FreeImage_Unload(gradients[k]); + gradients[k] = NULL; + } + } + return FALSE; + } +} + +/** +Compute the gradient attenuation function PHI(x, y) +@param gradients Gradient pyramid (nlevels levels) +@param avgGrad Average gradient on each level (array of size nlevels) +@param nlevels Number of levels +@param alpha Parameter alpha in the paper +@param beta Parameter beta in the paper +@return Returns the attenuation matrix Phi if successful, returns NULL otherwise +*/ +static FIBITMAP* PhiMatrix(FIBITMAP **gradients, float *avgGrad, int nlevels, float alpha, float beta) { + float *src_pixel, *dst_pixel; + FIBITMAP **phi = NULL; + + try { + phi = (FIBITMAP**)malloc(nlevels * sizeof(FIBITMAP*)); + if (!phi) throw(1); + memset(phi, 0, nlevels * sizeof(FIBITMAP*)); + + for(int k = nlevels-1; k >= 0; k--) { + // compute phi(k) + + FIBITMAP *Gk = gradients[k]; + + const unsigned width = FreeImage_GetWidth(Gk); + const unsigned height = FreeImage_GetHeight(Gk); + const unsigned pitch = FreeImage_GetPitch(Gk) / sizeof(float); + + // parameter alpha is 0.1 times the average gradient magnitude + // also, note the factor of 2**k in the denominator; + // that is there to correct for the fact that an average gradient avgGrad(H) over 2**k pixels + // in the original image will appear as a gradient grad(Hk) = 2**k*avgGrad(H) over a single pixel in Hk. + float ALPHA = alpha * avgGrad[k] * (float)((int)1 << k); + if(ALPHA == 0) ALPHA = EPSILON; + + phi[k] = FreeImage_AllocateT(FIT_FLOAT, width, height); + if (!phi[k]) throw(1); + + src_pixel = (float*)FreeImage_GetBits(Gk); + dst_pixel = (float*)FreeImage_GetBits(phi[k]); + for(unsigned y = 0; y < height; y++) { + for(unsigned x = 0; x < width; x++) { + // compute (alpha / grad) * (grad / alpha) ** beta + const float v = src_pixel[x] / ALPHA; + const float value = (float)pow((float)v, (float)(beta-1)); + dst_pixel[x] = (value > 1) ? 1 : value; + } + // next line + src_pixel += pitch; + dst_pixel += pitch; + } + + if(k < nlevels-1) { + // compute PHI(k) = L( PHI(k+1)) * phi(k) + FIBITMAP *L = FreeImage_Rescale(phi[k+1], width, height, FILTER_BILINEAR); + if (!L) throw(1); + + src_pixel = (float*)FreeImage_GetBits(L); + dst_pixel = (float*)FreeImage_GetBits(phi[k]); + for(unsigned y = 0; y < height; y++) { + for(unsigned x = 0; x < width; x++) { + dst_pixel[x] *= src_pixel[x]; + } + // next line + src_pixel += pitch; + dst_pixel += pitch; + } + + FreeImage_Unload(L); + + // PHI(k+1) is no longer needed + FreeImage_Unload(phi[k+1]); + phi[k+1] = NULL; + } + + // next level + } + + // get the final result and return + FIBITMAP *dst = phi[0]; + + free(phi); + + return dst; + + } catch(int) { + if(phi) { + for(int k = nlevels-1; k >= 0; k--) { + if(phi[k]) FreeImage_Unload(phi[k]); + } + free(phi); + } + return NULL; + } +} + +/** +Compute gradients in x and y directions, attenuate them with the attenuation matrix, +then compute the divergence div G from the attenuated gradient. +@param H Normalized luminance +@param PHI Attenuation matrix +@return Returns the divergence matrix if successful, returns NULL otherwise +*/ +static FIBITMAP* Divergence(FIBITMAP *H, FIBITMAP *PHI) { + FIBITMAP *Gx = NULL, *Gy = NULL, *divG = NULL; + float *phi, *h, *gx, *gy, *divg; + + try { + const FREE_IMAGE_TYPE image_type = FreeImage_GetImageType(H); + if(image_type != FIT_FLOAT) throw(1); + + const unsigned width = FreeImage_GetWidth(H); + const unsigned height = FreeImage_GetHeight(H); + + Gx = FreeImage_AllocateT(image_type, width, height); + if (!Gx) throw(1); + Gy = FreeImage_AllocateT(image_type, width, height); + if (!Gy) throw(1); + + const unsigned pitch = FreeImage_GetPitch(H) / sizeof(float); + + // perform gradient attenuation + + phi = (float*)FreeImage_GetBits(PHI); + h = (float*)FreeImage_GetBits(H); + gx = (float*)FreeImage_GetBits(Gx); + gy = (float*)FreeImage_GetBits(Gy); + + for(unsigned y = 0; y < height; y++) { + const unsigned s = (y+1 == height ? y : y+1); + for(unsigned x = 0; x < width; x++) { + const unsigned e = (x+1 == width ? x : x+1); + // forward difference + const unsigned index = y*pitch + x; + const float phi_xy = phi[index]; + const float h_xy = h[index]; + gx[x] = (h[y*pitch+e] - h_xy) * phi_xy; // [H(x+1, y) - H(x, y)] * PHI(x, y) + gy[x] = (h[s*pitch+x] - h_xy) * phi_xy; // [H(x, y+1) - H(x, y)] * PHI(x, y) + } + // next line + gx += pitch; + gy += pitch; + } + + // calculate the divergence + + divG = FreeImage_AllocateT(image_type, width, height); + if (!divG) throw(1); + + gx = (float*)FreeImage_GetBits(Gx); + gy = (float*)FreeImage_GetBits(Gy); + divg = (float*)FreeImage_GetBits(divG); + + for(unsigned y0 = 0; y0 < height; y0++) { + for(unsigned x = 0; x < width; x++) { + // backward difference approximation + // divG = Gx(x, y) - Gx(x-1, y) + Gy(x, y) - Gy(x, y-1) + const unsigned index = y0*pitch + x; + divg[index] = gx[index] + gy[index]; + if(x > 0) divg[index] -= gx[index-1]; + if(y0 > 0) divg[index] -= gy[index-pitch]; + } + } + + // no longer needed ... + FreeImage_Unload(Gx); + FreeImage_Unload(Gy); + + // return the divergence + return divG; + + } catch(int) { + if(Gx) FreeImage_Unload(Gx); + if(Gy) FreeImage_Unload(Gy); + if(divG) FreeImage_Unload(divG); + return NULL; + } +} + +/** +Given the luminance channel, find max & min luminance values, +normalize to range 0..100 and take the logarithm. +@param Y Image luminance +@return Returns the normalized luminance H if successful, returns NULL otherwise +*/ +static FIBITMAP* LogLuminance(FIBITMAP *Y) { + FIBITMAP *H = NULL; + + try { + // get the luminance channel + FIBITMAP *H = FreeImage_Clone(Y); + if (!H) throw(1); + + const unsigned width = FreeImage_GetWidth(H); + const unsigned height = FreeImage_GetHeight(H); + const unsigned pitch = FreeImage_GetPitch(H); + + // find max & min luminance values + float maxLum = -1e20F, minLum = 1e20F; + + BYTE *bits = (BYTE*)FreeImage_GetBits(H); + for(unsigned y = 0; y < height; y++) { + const float *pixel = (float*)bits; + for(unsigned x = 0; x < width; x++) { + const float value = pixel[x]; + maxLum = (maxLum < value) ? value : maxLum; // max Luminance in the scene + minLum = (minLum < value) ? minLum : value; // min Luminance in the scene + } + // next line + bits += pitch; + } + if(maxLum == minLum) throw(1); + + // normalize to range 0..100 and take the logarithm + const float scale = 100.F / (maxLum - minLum); + bits = (BYTE*)FreeImage_GetBits(H); + for(unsigned y0 = 0; y0 < height; y0++) { + float *pixel = (float*)bits; + for(unsigned x = 0; x < width; x++) { + const float value = (pixel[x] - minLum) * scale; + pixel[x] = log(value + EPSILON); + } + // next line + bits += pitch; + } + + return H; + + } catch(int) { + if(H) FreeImage_Unload(H); + return NULL; + } +} + +/** +Given a normalized luminance, perform exponentiation and recover the log compressed image +@param Y Input/Output luminance image +*/ +static void ExpLuminance(FIBITMAP *Y) { + const unsigned width = FreeImage_GetWidth(Y); + const unsigned height = FreeImage_GetHeight(Y); + const unsigned pitch = FreeImage_GetPitch(Y); + + BYTE *bits = (BYTE*)FreeImage_GetBits(Y); + for(unsigned y = 0; y < height; y++) { + float *pixel = (float*)bits; + for(unsigned x = 0; x < width; x++) { + pixel[x] = exp(pixel[x]) - EPSILON; + } + bits += pitch; + } +} + +// -------------------------------------------------------------------------- + +/** +Gradient Domain HDR tone mapping operator +@param Y Image luminance values +@param alpha Parameter alpha of the paper (suggested value is 0.1) +@param beta Parameter beta of the paper (suggested value is between 0.8 and 0.9) +@return returns the tone mapped luminance +*/ +static FIBITMAP* tmoFattal02(FIBITMAP *Y, float alpha, float beta) { + const unsigned MIN_PYRAMID_SIZE = 32; // minimun size (width or height) of the coarsest level of the pyramid + + FIBITMAP *H = NULL; + FIBITMAP **pyramid = NULL; + FIBITMAP **gradients = NULL; + FIBITMAP *phy = NULL; + FIBITMAP *divG = NULL; + FIBITMAP *U = NULL; + float *avgGrad = NULL; + + int k; + int nlevels = 0; + + try { + // get the normalized luminance + FIBITMAP *H = LogLuminance(Y); + if (!H) throw(1); + + // get the number of levels for the pyramid + const unsigned width = FreeImage_GetWidth(H); + const unsigned height = FreeImage_GetHeight(H); + unsigned minsize = MIN(width, height); + while(minsize >= MIN_PYRAMID_SIZE) { + nlevels++; + minsize /= 2; + } + + // create the Gaussian pyramid + pyramid = (FIBITMAP**)malloc(nlevels * sizeof(FIBITMAP*)); + if (!pyramid) throw(1); + memset(pyramid, 0, nlevels * sizeof(FIBITMAP*)); + + if (!GaussianPyramid(H, pyramid, nlevels)) throw(1); + + // calculate gradient magnitude and its average value on each pyramid level + gradients = (FIBITMAP**)malloc(nlevels * sizeof(FIBITMAP*)); + if (!gradients) throw(1); + memset(gradients, 0, nlevels * sizeof(FIBITMAP*)); + avgGrad = (float*)malloc(nlevels * sizeof(float)); + if (!avgGrad) throw(1); + + if (!GradientPyramid(pyramid, nlevels, gradients, avgGrad)) throw(1); + + // free the Gaussian pyramid + for(k = 0; k < nlevels; k++) { + if(pyramid[k]) FreeImage_Unload(pyramid[k]); + } + free(pyramid); pyramid = NULL; + + // compute the gradient attenuation function PHI(x, y) + phy = PhiMatrix(gradients, avgGrad, nlevels, alpha, beta); + if (!phy) throw(1); + + // free the gradient pyramid + for(k = 0; k < nlevels; k++) { + if(gradients[k]) FreeImage_Unload(gradients[k]); + } + free(gradients); gradients = NULL; + free(avgGrad); avgGrad = NULL; + + // compute gradients in x and y directions, attenuate them with the attenuation matrix, + // then compute the divergence div G from the attenuated gradient. + divG = Divergence(H, phy); + if (!divG) throw(1); + + // H & phy no longer needed + FreeImage_Unload(H); H = NULL; + FreeImage_Unload(phy); phy = NULL; + + // solve the PDE (Poisson equation) using a multigrid solver and 3 cycles + FIBITMAP *U = FreeImage_MultigridPoissonSolver(divG, 3); + if (!U) throw(1); + + FreeImage_Unload(divG); + + // perform exponentiation and recover the log compressed image + ExpLuminance(U); + + return U; + + } catch(int) { + if(H) FreeImage_Unload(H); + if(pyramid) { + for(int i = 0; i < nlevels; i++) { + if(pyramid[i]) FreeImage_Unload(pyramid[i]); + } + free(pyramid); + } + if(gradients) { + for(int i = 0; i < nlevels; i++) { + if(gradients[i]) FreeImage_Unload(gradients[i]); + } + free(gradients); + } + if(avgGrad) free(avgGrad); + if(phy) FreeImage_Unload(phy); + if(divG) FreeImage_Unload(divG); + if(U) FreeImage_Unload(U); + + return NULL; + } +} + +// ---------------------------------------------------------- +// Main algorithm +// ---------------------------------------------------------- + +/** +Apply the Gradient Domain High Dynamic Range Compression to a RGBF image and convert to 24-bit RGB +@param dib Input RGBF / RGB16 image +@param color_saturation Color saturation (s parameter in the paper) in [0.4..0.6] +@param attenuation Atenuation factor (beta parameter in the paper) in [0.8..0.9] +@return Returns a 24-bit RGB image if successful, returns NULL otherwise +*/ +FIBITMAP* DLL_CALLCONV +FreeImage_TmoFattal02(FIBITMAP *dib, double color_saturation, double attenuation) { + const float alpha = 0.1F; // parameter alpha = 0.1 + const float beta = (float)MAX(0.8, MIN(0.9, attenuation)); // parameter beta = [0.8..0.9] + const float s = (float)MAX(0.4, MIN(0.6, color_saturation));// exponent s controls color saturation = [0.4..0.6] + + FIBITMAP *src = NULL; + FIBITMAP *Yin = NULL; + FIBITMAP *Yout = NULL; + FIBITMAP *dst = NULL; + + if (!FreeImage_HasPixels(dib)) return NULL; + + try { + + // convert to RGBF + src = FreeImage_ConvertToRGBF(dib); + if (!src) throw(1); + + // get the luminance channel + Yin = ConvertRGBFToY(src); + if (!Yin) throw(1); + + // perform the tone mapping + Yout = tmoFattal02(Yin, alpha, beta); + if (!Yout) throw(1); + + // clip low and high values and normalize to [0..1] + //NormalizeY(Yout, 0.001F, 0.995F); + NormalizeY(Yout, 0, 1); + + // compress the dynamic range + + const unsigned width = FreeImage_GetWidth(src); + const unsigned height = FreeImage_GetHeight(src); + + const unsigned rgb_pitch = FreeImage_GetPitch(src); + const unsigned y_pitch = FreeImage_GetPitch(Yin); + + BYTE *bits = (BYTE*)FreeImage_GetBits(src); + BYTE *bits_yin = (BYTE*)FreeImage_GetBits(Yin); + BYTE *bits_yout = (BYTE*)FreeImage_GetBits(Yout); + + for(unsigned y = 0; y < height; y++) { + float *Lin = (float*)bits_yin; + float *Lout = (float*)bits_yout; + float *color = (float*)bits; + for(unsigned x = 0; x < width; x++) { + for(unsigned c = 0; c < 3; c++) { + *color = (Lin[x] > 0) ? pow(*color/Lin[x], s) * Lout[x] : 0; + color++; + } + } + bits += rgb_pitch; + bits_yin += y_pitch; + bits_yout += y_pitch; + } + + // not needed anymore + FreeImage_Unload(Yin); Yin = NULL; + FreeImage_Unload(Yout); Yout = NULL; + + // clamp image highest values to display white, then convert to 24-bit RGB + dst = ClampConvertRGBFTo24(src); + + // clean-up and return + FreeImage_Unload(src); src = NULL; + + // copy metadata from src to dst + FreeImage_CloneMetadata(dst, dib); + + return dst; + + } catch(int) { + if(src) FreeImage_Unload(src); + if(Yin) FreeImage_Unload(Yin); + if(Yout) FreeImage_Unload(Yout); + return NULL; + } +} -- cgit v1.2.3