pfstools-2.2.0.tgz: .../reinhard02/tmo_reinhard02.cpp

"Fossies" - the Fresh Open Source Software Archive

Member "pfstools-2.2.0/src/tmo/reinhard02/tmo_reinhard02.cpp" (12 Aug 2021, 13908 Bytes) of package /linux/privat/pfstools-2.2.0.tgz:

As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "tmo_reinhard02.cpp" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report:

2.1.0_vs_2.2.0.

1 /** 2 * @file tmo_reinhard02.cpp 3 * @brief Tone map luminance channel using Reinhard02 model 4 * 5 * Implementation courtesy of Erik Reinhard. 6 * 7 * Original source code note: 8 * Tonemap.c University of Utah / Erik Reinhard / October 2001 9 * 10 * File taken from: 11 * http://www.cs.utah.edu/~reinhard/cdrom/source.html 12 * 13 * Port to PFS tools library by Grzegorz Krawczyk <krawczyk@mpi-sb.mpg.de> 14 * 15 * $Id: tmo_reinhard02.cpp,v 1.4 2009/04/15 11:49:32 julians37 Exp $ 16 */ 17 18 #define _USE_MATH_DEFINES 19 #include <cmath> 20 21 #include <config.h> 22 23 #include <stdlib.h> 24 #include <stdio.h> 25 #include <math.h> 26 27 #include "pfstmo.h" 28 29 #ifdef HAVE_ZFFT 30 #include <fft.h> 31 #else 32 // if no zfft library compile approximate sollution 33 #define APPROXIMATE 34 #endif 35 36 //--- from defines.h 37 typedef struct { 38 int xmax, ymax; /* image dimensions */ 39 } CVTS; 40 41 typedef double COLOR[3]; /* red, green, blue (or X,Y,Z) */ 42 //--- end of defines.h 43 44 45 static int width, height, scale; 46 static COLOR **image; 47 static double ***convolved_image; 48 static double sigma_0, sigma_1; 49 double **luminance; 50 51 static double k = 1. / (2. * 1.4142136); 52 static double key = 0.18; 53 static double threshold = 0.05; 54 static double phi = 8.; 55 static double white = 1e20; 56 static int scale_low = 1; 57 static int scale_high = 43; // 1.6^8 = 43 58 static int range = 8; 59 static int use_scales = 0; 60 static int use_border = 0; 61 static CVTS cvts = {0, 0}; 62 63 static bool temporal_coherent; 64 65 #ifdef APPROXIMATE 66 extern int PyramidHeight; // set by build_pyramid, defines actual pyramid size 67 extern double V1 (int x, int y, int level); 68 extern void build_pyramid (double **luminance, int image_width, int image_height); 69 extern void clean_pyramid (); 70 #else 71 72 static zomplex **filter_fft; 73 static zomplex *image_fft, *coeff; 74 75 #define V1(x,y,i) (convolved_image[i][x][y]) 76 77 #endif 78 79 #define SIGMA_I(i) (sigma_0 + ((double)i/(double)range)*(sigma_1 - sigma_0)) 80 #define S_I(i) (exp (SIGMA_I(i))) 81 #define V2(x,y,i) (V1(x,y,i+1)) 82 #define ACTIVITY(x,y,i) ((V1(x,y,i) - V2(x,y,i)) / (((key * pow (2., phi))/(S_I(i)*S_I(i))) + V1(x,y,i))) 83 84 85 86 /** 87 * Used to achieve temporal coherence 88 */ 89 template<class T> 90 class TemporalSmoothVariable 91 { 92 // const int hist_length = 100; 93 T value; 94 95 T getThreshold( T luminance ) 96 { 97 return 0.01 * luminance; 98 } 99 100 public: 101 TemporalSmoothVariable() : value( -1 ) 102 { 103 } 104 105 void set( T new_value ) 106 { 107 if( value == -1 ) 108 value = new_value; 109 else { 110 T delta = new_value - value; 111 const T threshold = getThreshold( (new_value + value)/2 ); 112 if( delta > threshold ) delta = threshold; 113 else if( delta < -threshold ) delta = -threshold; 114 value += delta; 115 } 116 } 117 118 T get() const 119 { 120 return value; 121 } 122 }; 123 124 125 static TemporalSmoothVariable<double> avg_luminance, max_luminance; 126 127 128 /* 129 * Kaiser-Bessel stuff 130 */ 131 132 static double alpha = 2.; /* Kaiser-Bessel window parameter */ 133 static double bbeta; /* Will contain bessel (PI * alpha) */ 134 135 /* 136 * Modified zeroeth order bessel function of the first kind 137 */ 138 139 double bessel (double x) 140 { 141 const double f = 1e-9; 142 int n = 1; 143 double s = 1.; 144 double d = 1.; 145 146 double t; 147 148 while (d > f * s) 149 { 150 t = x / (2. * n); 151 n++; 152 d *= t * t; 153 s += d; 154 } 155 return s; 156 } 157 158 /* 159 * Initialiser for Kaiser-Bessel computations 160 */ 161 162 void compute_bessel () 163 { 164 bbeta = bessel (M_PI * alpha); 165 } 166 167 /* 168 * Kaiser-Bessel function with window length M and parameter beta = 2. 169 * Window length M = min (width, height) / 2 170 */ 171 172 double kaiserbessel (double x, double y, double M) 173 { 174 double d = 1. - ((x*x + y*y) / (M * M)); 175 if (d <= 0.) 176 return 0.; 177 return bessel (M_PI * alpha * sqrt (d)) / bbeta; 178 } 179 180 181 /* 182 * FFT functions 183 */ 184 #ifdef HAVE_ZFFT 185 186 void initialise_fft (int width, int height) 187 { 188 coeff = zfft2di (width, height, NULL); 189 } 190 191 void compute_fft (zomplex *array, int width, int height) 192 { 193 zfft2d (-1, width, height, array, width, coeff); 194 } 195 196 void compute_inverse_fft (zomplex *array, int width, int height) 197 { 198 zfft2d (1, width, height, array, width, coeff); 199 } 200 201 202 // Compute Gaussian blurred images 203 void gaussian_filter (zomplex *filter, double scale, double k ) 204 { 205 int x, y; 206 double x1, y1, s; 207 double a = 1. / (k * scale); 208 double c = 1. / 4.; 209 210 for (y = 0; y < cvts.ymax; y++) 211 { 212 y1 = (y >= cvts.ymax / 2) ? y - cvts.ymax : y; 213 s = erf (a * (y1 - .5)) - erf (a * (y1 + .5)); 214 for (x = 0; x < cvts.xmax; x++) 215 { 216 x1 = (x >= cvts.xmax / 2) ? x - cvts.xmax : x; 217 filter[y*cvts.xmax + x].re = s * (erf (a * (x1 - .5)) - erf (a * (x1 + .5))) * c; 218 filter[y*cvts.xmax + x].im = 0.; 219 } 220 } 221 } 222 223 void build_gaussian_fft () 224 { 225 int i; 226 double length = cvts.xmax * cvts.ymax; 227 double fft_scale = 1. / sqrt (length); 228 filter_fft = (zomplex**) calloc (range, sizeof (zomplex*)); 229 230 for (scale = 0; scale < range; scale++) 231 { 232 fprintf (stderr, "Computing FFT of Gaussian at scale %i (size %i x %i)%c", 233 scale, cvts.xmax, cvts.ymax, (char)13); 234 filter_fft[scale] = (zomplex*) calloc (length, sizeof (zomplex)); 235 gaussian_filter (filter_fft[scale], S_I(scale), k); 236 compute_fft (filter_fft[scale], cvts.xmax, cvts.ymax); 237 for (i = 0; i < length; i++) 238 { 239 filter_fft[scale][i].re *= fft_scale; 240 filter_fft[scale][i].im *= fft_scale; 241 } 242 } 243 fprintf (stderr, "\n"); 244 } 245 246 void build_image_fft () 247 { 248 int i, x, y; 249 double length = cvts.xmax * cvts.ymax; 250 double fft_scale = 1. / sqrt (length); 251 252 fprintf (stderr, "Computing image FFT\n"); 253 image_fft = (zomplex*) calloc (length, sizeof (zomplex)); 254 255 for (y = 0; y < cvts.ymax; y++) 256 for (x = 0; x < cvts.xmax; x++) 257 image_fft[y*cvts.xmax + x].re = luminance[y][x]; 258 259 compute_fft (image_fft, cvts.xmax, cvts.ymax); 260 for (i = 0; i < length; i++) 261 { 262 image_fft[i].re *= fft_scale; 263 image_fft[i].im *= fft_scale; 264 } 265 } 266 267 void convolve_filter (int scale, zomplex *convolution_fft) 268 { 269 int i, x, y; 270 271 for (i = 0; i < cvts.xmax * cvts.ymax; i++) 272 { 273 convolution_fft[i].re = image_fft[i].re * filter_fft[scale][i].re - 274 image_fft[i].im * filter_fft[scale][i].im; 275 convolution_fft[i].im = image_fft[i].re * filter_fft[scale][i].im + 276 image_fft[i].im * filter_fft[scale][i].re; 277 } 278 compute_inverse_fft (convolution_fft, cvts.xmax, cvts.ymax); 279 i = 0; 280 for (y = 0; y < cvts.ymax; y++) 281 for (x = 0; x < cvts.xmax; x++) 282 convolved_image[scale][x][y] = convolution_fft[i++].re; 283 } 284 285 void compute_fourier_convolution () 286 { 287 int x; 288 zomplex *convolution_fft; 289 290 initialise_fft (cvts.xmax, cvts.ymax); 291 build_image_fft (); 292 build_gaussian_fft (); 293 convolved_image = (double ***) malloc (range * sizeof (double **)); 294 295 convolution_fft = (zomplex *) calloc (cvts.xmax * cvts.ymax, sizeof (zomplex)); 296 for (scale = 0; scale < range; scale++) 297 { 298 fprintf (stderr, "Computing convolved image at scale %i%c", scale, (char)13); 299 convolved_image[scale] = (double **) malloc (cvts.xmax * sizeof (double *)); 300 for (x = 0; x < cvts.xmax; x++) 301 convolved_image[scale][x] = (double *) malloc (cvts.ymax * sizeof (double)); 302 convolve_filter (scale, convolution_fft); 303 free (filter_fft[scale]); 304 } 305 fprintf (stderr, "\n"); 306 free (convolution_fft); 307 free (image_fft); 308 } 309 310 #endif // #ifdef HAVE_ZFFT 311 312 313 314 /* 315 * Tonemapping routines 316 */ 317 318 double get_maxvalue () 319 { 320 double max = 0.; 321 int x, y; 322 323 for (y = 0; y < cvts.ymax; y++) 324 for (x = 0; x < cvts.xmax; x++) 325 max = (max < image[y][x][0]) ? image[y][x][0] : max; 326 return max; 327 } 328 329 void tonemap_image () 330 { 331 double Lmax2; 332 int x, y; 333 int scale, prefscale; 334 335 if (white < 1e20) 336 Lmax2 = white * white; 337 else 338 { 339 if( temporal_coherent ) { 340 max_luminance.set( get_maxvalue() ); 341 Lmax2 = max_luminance.get(); 342 } else Lmax2 = get_maxvalue(); 343 Lmax2 *= Lmax2; 344 } 345 346 for (y = 0; y < cvts.ymax; y++) 347 for (x = 0; x < cvts.xmax; x++) 348 { 349 if (use_scales) 350 { 351 prefscale = range - 1; 352 for (scale = 0; scale < range - 1; scale++) 353 if ( scale >= PyramidHeight || fabs (ACTIVITY(x,y,scale)) > threshold) 354 { 355 prefscale = scale; 356 break; 357 } 358 image[y][x][0] /= 1. + V1(x,y,prefscale); 359 } 360 else 361 image[y][x][0] = image[y][x][0] * (1. + (image[y][x][0] / Lmax2)) / (1. + image[y][x][0]); 362 // image[y][x][0] /= (1. + image[y][x][0]); 363 } 364 } 365 366 /* 367 * Miscellaneous functions 368 */ 369 370 void clamp_image () 371 { 372 int x, y; 373 374 for (y = 0; y < cvts.ymax; y++) 375 for (x = 0; x < cvts.xmax; x++) 376 { 377 image[y][x][0] = (image[y][x][0] > 1.) ? 1. : image[y][x][0]; 378 image[y][x][1] = (image[y][x][1] > 1.) ? 1. : image[y][x][1]; 379 image[y][x][2] = (image[y][x][2] > 1.) ? 1. : image[y][x][2]; 380 } 381 } 382 383 double log_average () 384 { 385 int x, y; 386 double sum = 0.; 387 388 for (y = 0; y < cvts.ymax; y++) 389 for (x = 0; x < cvts.xmax; x++) 390 sum += log (0.00001 + luminance[y][x]); 391 return exp (sum / (double)(cvts.xmax * cvts.ymax)); 392 } 393 394 void scale_to_midtone () 395 { 396 int x, y, u, v, d; 397 double factor; 398 double scale_factor; 399 double low_tone = key / 3.; 400 int border_size = (cvts.xmax < cvts.ymax) ? int(cvts.xmax / 5.) : int(cvts.ymax / 5.); 401 int hw = cvts.xmax >> 1; 402 int hh = cvts.ymax >> 1; 403 404 double avg; 405 if( temporal_coherent ) { 406 avg_luminance.set( log_average() ); 407 avg = avg_luminance.get(); 408 } else avg = log_average(); 409 410 scale_factor = 1.0 / avg; 411 for (y = 0; y < cvts.ymax; y++) 412 for (x = 0; x < cvts.xmax; x++) 413 { 414 if (use_border) 415 { 416 u = (x > hw) ? cvts.xmax - x : x; 417 v = (y > hh) ? cvts.ymax - y : y; 418 d = (u < v) ? u : v; 419 factor = (d < border_size) ? (key - low_tone) * 420 kaiserbessel (border_size - d, 0, border_size) + 421 low_tone : key; 422 } 423 else 424 factor = key; 425 image[y][x][0] *= scale_factor * factor; 426 luminance[y][x] *= scale_factor * factor; 427 } 428 } 429 430 void copy_luminance () 431 { 432 int x, y; 433 434 for (x = 0; x < cvts.xmax; x++) 435 for (y = 0; y < cvts.ymax; y++) 436 luminance[y][x] = image[y][x][0]; 437 } 438 439 /* 440 * Memory allocation 441 */ 442 void allocate_memory () 443 { 444 int y; 445 446 luminance = (double **) malloc (cvts.ymax * sizeof (double *)); 447 image = (COLOR **) malloc (cvts.ymax * sizeof (COLOR *)); 448 for (y = 0; y < cvts.ymax; y++) 449 { 450 luminance[y] = (double *) malloc (cvts.xmax * sizeof (double)); 451 image[y] = (COLOR *) malloc (cvts.xmax * sizeof (COLOR)); 452 } 453 } 454 455 void deallocate_memory () 456 { 457 int y; 458 459 for (y = 0; y < cvts.ymax; y++) 460 { 461 free(luminance[y]); 462 free(image[y]); 463 } 464 free( luminance ); 465 free( image ); 466 } 467 468 469 void dynamic_range () 470 { 471 int x, y; 472 double minval = 1e20; 473 double maxval = -1e20; 474 475 for (x = 0; x < cvts.xmax; x++) 476 for (y = 0; y < cvts.ymax; y++) 477 { 478 if ((luminance[y][x] < minval) && 479 (luminance[y][x] > 0.0)) 480 minval = luminance[y][x]; 481 if (luminance[y][x] > maxval) 482 maxval = luminance[y][x]; 483 } 484 fprintf (stderr, "\tRange of values = %9.8f - %9.8f\n", minval, maxval); 485 fprintf (stderr, "\tDynamic range = %i:1\n", (int)(maxval/minval)); 486 } 487 488 void print_parameter_settings () 489 { 490 fprintf (stderr, "\tImage size = %i %i\n", cvts.xmax, cvts.ymax); 491 fprintf (stderr, "\tLowest scale = %i pixels\t\t(-low <integer>)\n", scale_low); 492 fprintf (stderr, "\tHighest scale = %i pixels\t\t(-high <integer>)\n", scale_high); 493 fprintf (stderr, "\tNumber of scales = %i\t\t\t(-num <integer>)\n", range); 494 fprintf (stderr, "\tScale spacing = %f\n", S_I(1) / S_I(0)); 495 fprintf (stderr, "\tKey value = %f\t\t(-key <float>)\n", key); 496 fprintf (stderr, "\tWhite value = %f\t\t(-white <float>)\n", white); 497 fprintf (stderr, "\tPhi = %f\t\t(-phi <float>)\n", phi); 498 fprintf (stderr, "\tThreshold = %f\t\t(-threshold <float>)\n", threshold); 499 dynamic_range (); 500 } 501 502 /* 503 * @brief Photographic tone-reproduction 504 * 505 * @param Y input luminance 506 * @param L output tonemapped intensities 507 * @param use_scales true: local version, false: global version of TMO 508 * @param key maps log average luminance to this value (default: 0.18) 509 * @param phi sharpening parameter (defaults to 1 - no sharpening) 510 * @param num number of scales to use in computation (default: 8) 511 * @param low size in pixels of smallest scale (should be kept at 1) 512 * @param high size in pixels of largest scale (default 1.6^8 = 43) 513 */ 514 void tmo_reinhard02( 515 unsigned int width, unsigned int height, 516 const float *nY, float *nL, 517 bool use_scales, float key, float phi, 518 int num, int low, int high, bool temporal_coherent ) 519 { 520 const pfstmo::Array2D* Y = new pfstmo::Array2D(width, height, const_cast<float*>(nY)); 521 pfstmo::Array2D* L = new pfstmo::Array2D(width, height, nL); 522 523 int x,y; 524 525 ::key = key; 526 ::phi = phi; 527 ::range = num; 528 ::scale_low = low; 529 ::scale_high = high; 530 ::use_scales = (use_scales) ? 1 : 0; 531 ::temporal_coherent = temporal_coherent; 532 533 cvts.xmax = Y->getCols(); 534 cvts.ymax = Y->getRows(); 535 536 sigma_0 = log (scale_low); 537 sigma_1 = log (scale_high); 538 539 compute_bessel(); 540 allocate_memory (); 541 542 // reading image 543 for( y=0 ; y<cvts.ymax ; y++ ) 544 for( x=0 ; x<cvts.xmax ; x++ ) 545 image[y][x][0] = (*Y)(x,y); 546 547 copy_luminance(); 548 scale_to_midtone(); 549 550 if( use_scales ) 551 { 552 #ifdef APPROXIMATE 553 build_pyramid(luminance, cvts.xmax, cvts.ymax); 554 #else 555 compute_fourier_convolution(); 556 #endif 557 } 558 559 tonemap_image(); 560 561 // saving image 562 for( y=0 ; y<cvts.ymax ; y++ ) 563 for( x=0 ; x<cvts.xmax ; x++ ) 564 (*L)(x,y) = image[y][x][0]; 565 566 // print_parameter_settings(); 567 568 deallocate_memory(); 569 clean_pyramid(); 570 571 delete L; 572 delete Y; 573 }