"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.
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 }