"Fossies" - the Fresh Open Source Software Archive 
"Fossies" - the Fresh Open Source Software Archive A hint: This file contains one or more very long lines, so maybe it is better readable using the pure text view mode that shows the contents as wrapped lines within the browser window.
1 /** 2 * @Brief Contrast mapping TMO 3 * 4 * From: 5 * 6 * Rafal Mantiuk, Karol Myszkowski, Hans-Peter Seidel. 7 * A Perceptual Framework for Contrast Processing of High Dynamic Range Images 8 * In: ACM Transactions on Applied Perception 3 (3), pp. 286-308, 2006 9 * http://www.mpi-inf.mpg.de/~mantiuk/contrast_domain/ 10 * 11 * This file is a part of PFSTMO package. 12 * ---------------------------------------------------------------------- 13 * Copyright (C) 2007 Grzegorz Krawczyk 14 * 15 * This program is free software; you can redistribute it and/or modify 16 * it under the terms of the GNU General Public License as published by 17 * the Free Software Foundation; either version 2 of the License, or 18 * (at your option) any later version. 19 * 20 * This program is distributed in the hope that it will be useful, 21 * but WITHOUT ANY WARRANTY; without even the implied warranty of 22 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 23 * GNU General Public License for more details. 24 * 25 * You should have received a copy of the GNU General Public License 26 * along with this program; if not, write to the Free Software 27 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 28 * ---------------------------------------------------------------------- 29 * 30 * @author Radoslaw Mantiuk, <radoslaw.mantiuk@gmail.com> 31 * @author Rafal Mantiuk, <mantiuk@gmail.com> 32 * Updated 2007/12/17 by Ed Brambley <E.J.Brambley@damtp.cam.ac.uk> 33 * (more information on the changes: 34 * http://www.damtp.cam.ac.uk/user/ejb48/hdr/index.html) 35 * Updated 2008/06/25 by Ed Brambley <E.J.Brambley@damtp.cam.ac.uk> 36 * bug fixes and openMP patches 37 * more on this: 38 * http://groups.google.com/group/pfstools/browse_thread/thread/de2378af98ec6185/0dee5304fc14e99d?hl=en#0dee5304fc14e99d 39 * Optimization improvements by Lebed Dmytry 40 * 41 * $Id: contrast_domain.cpp,v 1.16 2012/04/22 13:14:14 rafm Exp $ 42 */ 43 44 #include <stdio.h> 45 #include <stdlib.h> 46 #include <string.h> 47 #include <math.h> 48 49 #include "contrast_domain.h" 50 51 #include "config.h" 52 53 typedef struct pyramid_s { 54 int rows; 55 int cols; 56 float* Gx; 57 float* Gy; 58 struct pyramid_s* next; 59 struct pyramid_s* prev; 60 } pyramid_t; 61 62 63 extern float xyz2rgbD65Mat[3][3]; 64 extern float rgb2xyzD65Mat[3][3]; 65 66 #define PYRAMID_MIN_PIXELS 3 67 #define LOOKUP_W_TO_R 107 68 69 static void contrast_equalization( pyramid_t *pp, const float contrastFactor ); 70 71 static void transform_to_luminance(pyramid_t* pyramid, float* const x, pfstmo_progress_callback progress_cb, const bool bcg); 72 static void matrix_add(const int n, const float* const a, float* const b); 73 static void matrix_subtract(const int n, const float* const a, float* const b); 74 static void matrix_copy(const int n, const float* const a, float* const b); 75 static void matrix_multiply_const(const int n, float* const a, const float val); 76 static void matrix_divide(const int n, const float* const a, float* const b); 77 static float* matrix_alloc(const int size); 78 static void matrix_free(float* m); 79 static float matrix_DotProduct(const int n, const float* const a, const float* const b); 80 static void matrix_zero(const int n, float* const m); 81 static void calculate_and_add_divergence(const int rows, const int cols, const float* const Gx, const float* const Gy, float* const divG); 82 static void pyramid_calculate_divergence(pyramid_t* pyramid); 83 static void pyramid_calculate_divergence_sum(pyramid_t* pyramid, float* divG_sum); 84 static void calculate_scale_factor(const int n, const float* const G, float* const C); 85 static void pyramid_calculate_scale_factor(pyramid_t* pyramid, pyramid_t* pC); 86 static void scale_gradient(const int n, float* const G, const float* const C); 87 static void pyramid_scale_gradient(pyramid_t* pyramid, pyramid_t* pC); 88 static void pyramid_free(pyramid_t* pyramid); 89 static pyramid_t* pyramid_allocate(const int cols, const int rows); 90 static void calculate_gradient(const int cols, const int rows, const float* const lum, float* const Gx, float* const Gy); 91 static void pyramid_calculate_gradient(pyramid_t* pyramid, float* lum); 92 static void solveX(const int n, const float* const b, float* const x); 93 static void multiplyA(pyramid_t* px, pyramid_t* pyramid, const float* const x, float* const divG_sum); 94 static void linbcg(pyramid_t* pyramid, pyramid_t* pC, const float* const b, float* const x, const int itmax, const float tol, pfstmo_progress_callback progress_cb); 95 static void lincg(pyramid_t* pyramid, pyramid_t* pC, const float* const b, float* const x, const int itmax, const float tol, pfstmo_progress_callback progress_cb); 96 static float lookup_table(const int n, const float* const in_tab, const float* const out_tab, const float val); 97 static void transform_to_R(const int n, float* const G); 98 static void pyramid_transform_to_R(pyramid_t* pyramid); 99 static void transform_to_G(const int n, float* const R); 100 static void pyramid_transform_to_G(pyramid_t* pyramid); 101 static void pyramid_gradient_multiply(pyramid_t* pyramid, const float val); 102 103 static void dump_matrix_to_file(const int width, const int height, const float* const m, const char * const file_name); 104 static void matrix_show(const char* const text, int rows, int cols, const float* const data); 105 static void pyramid_show(pyramid_t* pyramid); 106 107 108 static float W_table[] = {0.000000,0.010000,0.021180,0.031830,0.042628,0.053819,0.065556,0.077960,0.091140,0.105203,0.120255,0.136410,0.153788,0.172518,0.192739,0.214605,0.238282,0.263952,0.291817,0.322099,0.355040,0.390911,0.430009,0.472663,0.519238,0.570138,0.625811,0.686754,0.753519,0.826720,0.907041,0.995242,1.092169,1.198767,1.316090,1.445315,1.587756,1.744884,1.918345,2.109983,2.321863,2.556306,2.815914,3.103613,3.422694,3.776862,4.170291,4.607686,5.094361,5.636316,6.240338,6.914106,7.666321,8.506849,9.446889,10.499164,11.678143,13.000302,14.484414,16.151900,18.027221,20.138345,22.517282,25.200713,28.230715,31.655611,35.530967,39.920749,44.898685,50.549857,56.972578,64.280589,72.605654,82.100619,92.943020,105.339358,119.530154,135.795960,154.464484,175.919088,200.608905,229.060934,261.894494,299.838552,343.752526,394.651294,453.735325,522.427053,602.414859,695.706358,804.693100,932.229271,1081.727632,1257.276717,1463.784297,1707.153398,1994.498731,2334.413424,2737.298517,3215.770944,3785.169959,4464.187290,5275.653272,6247.520102,7414.094945,8817.590551,10510.080619}; 109 static float R_table[] = {0.000000,0.009434,0.018868,0.028302,0.037736,0.047170,0.056604,0.066038,0.075472,0.084906,0.094340,0.103774,0.113208,0.122642,0.132075,0.141509,0.150943,0.160377,0.169811,0.179245,0.188679,0.198113,0.207547,0.216981,0.226415,0.235849,0.245283,0.254717,0.264151,0.273585,0.283019,0.292453,0.301887,0.311321,0.320755,0.330189,0.339623,0.349057,0.358491,0.367925,0.377358,0.386792,0.396226,0.405660,0.415094,0.424528,0.433962,0.443396,0.452830,0.462264,0.471698,0.481132,0.490566,0.500000,0.509434,0.518868,0.528302,0.537736,0.547170,0.556604,0.566038,0.575472,0.584906,0.594340,0.603774,0.613208,0.622642,0.632075,0.641509,0.650943,0.660377,0.669811,0.679245,0.688679,0.698113,0.707547,0.716981,0.726415,0.735849,0.745283,0.754717,0.764151,0.773585,0.783019,0.792453,0.801887,0.811321,0.820755,0.830189,0.839623,0.849057,0.858491,0.867925,0.877358,0.886792,0.896226,0.905660,0.915094,0.924528,0.933962,0.943396,0.952830,0.962264,0.971698,0.981132,0.990566,1.000000}; 110 111 112 // display matrix in the console (debugging) 113 static void matrix_show(const char* const text, int cols, int rows, const float* const data) 114 { 115 const int _cols = cols; 116 117 if(rows > 8) 118 rows = 8; 119 if(cols > 8) 120 cols = 8; 121 122 printf("\n%s\n", text); 123 for(int ky=0; ky<rows; ky++){ 124 for(int kx=0; kx<cols; kx++){ 125 printf("%.06f ", data[kx + ky*_cols]); 126 } 127 printf("\n"); 128 } 129 } 130 131 132 // display pyramid in the console (debugging) 133 static void pyramid_show(pyramid_t* pyramid) 134 { 135 char ss[30]; 136 137 while (pyramid->next != NULL) 138 pyramid = pyramid->next; 139 140 while (pyramid != NULL) 141 { 142 printf("\n----- pyramid_t level %d,%d\n", pyramid->cols, pyramid->rows); 143 144 sprintf(ss, "Gx %p ", pyramid->Gx); 145 if(pyramid->Gx != NULL) 146 matrix_show(ss,pyramid->cols, pyramid->rows, pyramid->Gx); 147 sprintf(ss, "Gy %p ", pyramid->Gy); 148 if(pyramid->Gy != NULL) 149 matrix_show(ss,pyramid->cols, pyramid->rows, pyramid->Gy); 150 151 pyramid = pyramid->prev; 152 } 153 } 154 155 156 157 static inline float max( float a, float b ) 158 { 159 return a > b ? a : b; 160 } 161 162 static inline float min( float a, float b ) 163 { 164 return a < b ? a : b; 165 } 166 167 static inline int imin(int a, int b) 168 { 169 return a < b ? a : b; 170 } 171 172 173 // upsample the matrix 174 // upsampled matrix is twice bigger in each direction than data[] 175 // res should be a pointer to allocated memory for bigger matrix 176 // cols and rows are the dimmensions of the output matrix 177 static void matrix_upsample(const int outCols, const int outRows, const float* const in, float* const out) 178 { 179 const int inRows = outRows/2; 180 const int inCols = outCols/2; 181 182 // Transpose of experimental downsampling matrix (theoretically the correct thing to do) 183 184 const float dx = (float)inCols / ((float)outCols); 185 const float dy = (float)inRows / ((float)outRows); 186 const float factor = 1.0f / (dx*dy); // This gives a genuine upsampling matrix, not the transpose of the downsampling matrix 187 // const float factor = 1.0f; // Theoretically, this should be the best. 188 189 #pragma omp parallel for schedule(static) 190 for (int y = 0; y < outRows; y++) 191 { 192 const float sy = y * dy; 193 const int iy1 = ( y * inRows) / outRows; 194 const int iy2 = imin(((y+1) * inRows) / outRows, inRows-1); 195 196 for (int x = 0; x < outCols; x++) 197 { 198 const float sx = x * dx; 199 const int ix1 = ( x * inCols) / outCols; 200 const int ix2 = imin(((x+1) * inCols) / outCols, inCols-1); 201 202 out[x + y*outCols] = ( 203 ((ix1+1) - sx)*((iy1+1 - sy)) * in[ix1 + iy1*inCols] + 204 ((ix1+1) - sx)*(sy+dy - (iy1+1)) * in[ix1 + iy2*inCols] + 205 (sx+dx - (ix1+1))*((iy1+1 - sy)) * in[ix2 + iy1*inCols] + 206 (sx+dx - (ix1+1))*(sy+dx - (iy1+1)) * in[ix2 + iy2*inCols])*factor; 207 } 208 } 209 } 210 211 212 // downsample the matrix 213 static void matrix_downsample(const int inCols, const int inRows, const float* const data, float* const res) 214 { 215 const int outRows = inRows / 2; 216 const int outCols = inCols / 2; 217 218 const float dx = (float)inCols / ((float)outCols); 219 const float dy = (float)inRows / ((float)outRows); 220 221 // New downsampling by Ed Brambley: 222 // Experimental downsampling that assumes pixels are square and 223 // integrates over each new pixel to find the average value of the 224 // underlying pixels. 225 // 226 // Consider the original pixels laid out, and the new (larger) 227 // pixels layed out over the top of them. Then the new value for 228 // the larger pixels is just the integral over that pixel of what 229 // shows through; i.e., the values of the pixels underneath 230 // multiplied by how much of that pixel is showing. 231 // 232 // (ix1, iy1) is the coordinate of the top left visible pixel. 233 // (ix2, iy2) is the coordinate of the bottom right visible pixel. 234 // (fx1, fy1) is the fraction of the top left pixel showing. 235 // (fx2, fy2) is the fraction of the bottom right pixel showing. 236 237 const float normalize = 1.0f/(dx*dy); 238 #pragma omp parallel for schedule(static) 239 for (int y = 0; y < outRows; y++) 240 { 241 const int iy1 = ( y * inRows) / outRows; 242 const int iy2 = ((y+1) * inRows) / outRows; 243 const float fy1 = (iy1+1) - y * dy; 244 const float fy2 = (y+1) * dy - iy2; 245 246 for (int x = 0; x < outCols; x++) 247 { 248 const int ix1 = ( x * inCols) / outCols; 249 const int ix2 = ((x+1) * inCols) / outCols; 250 const float fx1 = (ix1+1) - x * dx; 251 const float fx2 = (x+1) * dx - ix2; 252 253 float pixVal = 0.0f; 254 float factorx, factory; 255 for (int i = iy1; i <= iy2 && i < inRows; i++) 256 { 257 if (i == iy1) 258 factory = fy1; // We're just getting the bottom edge of this pixel 259 else if (i == iy2) 260 factory = fy2; // We're just gettting the top edge of this pixel 261 else 262 factory = 1.0f; // We've got the full height of this pixel 263 for (int j = ix1; j <= ix2 && j < inCols; j++) 264 { 265 if (j == ix1) 266 factorx = fx1; // We've just got the right edge of this pixel 267 else if (j == ix2) 268 factorx = fx2; // We've just got the left edge of this pixel 269 else 270 factorx = 1.0f; // We've got the full width of this pixel 271 272 pixVal += data[j + i*inCols] * factorx * factory; 273 } 274 } 275 276 res[x + y * outCols] = pixVal * normalize; // Normalize by the area of the new pixel 277 } 278 } 279 } 280 281 // return = a + b 282 static inline void matrix_add(const int n, const float* const a, float* const b) 283 { 284 #pragma omp parallel for schedule(static) 285 for(int i=0; i<n; i++) 286 b[i] += a[i]; 287 } 288 289 // return = a - b 290 static inline void matrix_subtract(const int n, const float* const a, float* const b) 291 { 292 #pragma omp parallel for schedule(static) 293 for(int i=0; i<n; i++) 294 b[i] = a[i] - b[i]; 295 } 296 297 // copy matix a to b, return = a 298 static inline void matrix_copy(const int n, const float* const a, float* const b) 299 { 300 memcpy(b, a, sizeof(float)*n); 301 } 302 303 // multiply matrix a by scalar val 304 static inline void matrix_multiply_const(const int n, float* const a, const float val) 305 { 306 #pragma omp parallel for schedule(static) 307 for(int i=0; i<n; i++) 308 a[i] *= val; 309 } 310 311 // b = a[i] / b[i] 312 static inline void matrix_divide(const int n, const float* const a, float* const b) 313 { 314 #pragma omp parallel for schedule(static) 315 for(int i=0; i<n; i++) 316 b[i] = a[i] / b[i]; 317 } 318 319 320 // alloc memory for the float table 321 static inline float* matrix_alloc(int size){ 322 323 float* m = (float*)malloc(sizeof(float)*size); 324 if(m == NULL){ 325 fprintf(stderr, "ERROR: malloc in matrix_alloc() (size:%d)", size); 326 exit(155); 327 } 328 329 return m; 330 } 331 332 // free memory for matrix 333 static inline void matrix_free(float* m){ 334 if(m != NULL) 335 free(m); 336 } 337 338 // multiply vector by vector (each vector should have one dimension equal to 1) 339 static inline float matrix_DotProduct(const int n, const float* const a, const float* const b){ 340 float val = 0; 341 342 #pragma omp parallel for reduction(+:val) schedule(static) 343 for(int j=0;j<n;j++) 344 val += a[j] * b[j]; 345 346 return val; 347 } 348 349 // set zeros for matrix elements 350 static inline void matrix_zero(int n, float* m) 351 { 352 //bzero(m, n*sizeof(float)); 353 memset(m, n, sizeof(float)); 354 } 355 356 // calculate divergence of two gradient maps (Gx and Gy) 357 // divG(x,y) = Gx(x,y) - Gx(x-1,y) + Gy(x,y) - Gy(x,y-1) 358 static inline void calculate_and_add_divergence(const int cols, const int rows, const float* const Gx, const float* const Gy, float* const divG) 359 { 360 #pragma omp parallel for schedule(static) 361 for(int ky=0; ky<rows; ky++) 362 for(int kx=0; kx<cols; kx++) 363 { 364 float divGx, divGy; 365 const int idx = kx + ky*cols; 366 367 if(kx == 0) 368 divGx = Gx[idx]; 369 else 370 divGx = Gx[idx] - Gx[idx-1]; 371 372 if(ky == 0) 373 divGy = Gy[idx]; 374 else 375 divGy = Gy[idx] - Gy[idx - cols]; 376 377 divG[idx] += divGx + divGy; 378 } 379 } 380 381 // calculate the sum of divergences for the all pyramid level 382 // the smaller divergence map is upsamled and added to the divergence map for the higher level of pyramid 383 // temp is a temporary matrix of size (cols, rows), assumed to already be allocated 384 static void pyramid_calculate_divergence_sum(pyramid_t* pyramid, float* divG_sum) 385 { 386 float* temp = matrix_alloc(pyramid->rows*pyramid->cols); 387 388 // Find the coarsest pyramid, and the number of pyramid levels 389 int levels = 1; 390 while (pyramid->next != NULL) 391 { 392 levels++; 393 pyramid = pyramid->next; 394 } 395 396 // For every level, we swap temp and divG_sum. So, if there are an odd number of levels... 397 if (levels % 2) 398 { 399 float* const dummy = divG_sum; 400 divG_sum = temp; 401 temp = dummy; 402 } 403 404 // Add them all together 405 while (pyramid != NULL) 406 { 407 // Upsample or zero as needed 408 if (pyramid->next != NULL) 409 matrix_upsample(pyramid->cols, pyramid->rows, divG_sum, temp); 410 else 411 matrix_zero(pyramid->rows * pyramid->cols, temp); 412 413 // Add in the (freshly calculated) divergences 414 calculate_and_add_divergence(pyramid->cols, pyramid->rows, pyramid->Gx, pyramid->Gy, temp); 415 416 // char name[256]; 417 // sprintf( name, "Up_%d.pfs", pyramid->cols ); 418 // dump_matrix_to_file( pyramid->cols, pyramid->rows, temp, name ); 419 420 // matrix_copy(pyramid->rows*pyramid->cols, temp, divG_sum); 421 422 // Rather than copying, just switch round the pointers: we know we get them the right way round at the end. 423 float* const dummy = divG_sum; 424 divG_sum = temp; 425 temp = dummy; 426 427 pyramid = pyramid->prev; 428 } 429 430 matrix_free(temp); 431 } 432 433 // calculate scale factors (Cx,Cy) for gradients (Gx,Gy) 434 // C is equal to EDGE_WEIGHT for gradients smaller than GFIXATE or 1.0 otherwise 435 static inline void calculate_scale_factor(const int n, const float* const G, float* const C) 436 { 437 float GFIXATE = 0.1f; 438 float EDGE_WEIGHT = 0.01f; 439 const float detectT = 0.001f; 440 const float a = 0.038737; 441 const float b = 0.537756; 442 443 #pragma omp parallel for schedule(static) 444 for(int i=0; i<n; i++) 445 { 446 #if 1 447 const float g = max( detectT, fabsf(G[i]) ); 448 C[i] = 1.0f / (a*powf(g,b)); 449 #else 450 if(fabsf(G[i]) < GFIXATE) 451 C[i] = 1.0f / EDGE_WEIGHT; 452 else 453 C[i] = 1.0f; 454 #endif 455 } 456 } 457 458 // calculate scale factor for the whole pyramid 459 static void pyramid_calculate_scale_factor(pyramid_t* pyramid, pyramid_t* pC) 460 { 461 while (pyramid != NULL) 462 { 463 const int size = pyramid->rows * pyramid->cols; 464 calculate_scale_factor(size, pyramid->Gx, pC->Gx); 465 calculate_scale_factor(size, pyramid->Gy, pC->Gy); 466 pyramid = pyramid->next; 467 pC = pC->next; 468 } 469 } 470 471 // Scale gradient (Gx and Gy) by C (Cx and Cy) 472 // G = G / C 473 static inline void scale_gradient(const int n, float* const G, const float* const C) 474 { 475 #pragma omp parallel for schedule(static) 476 for(int i=0; i<n; i++) 477 G[i] *= C[i]; 478 } 479 480 // scale gradients for the whole one pyramid with the use of (Cx,Cy) from the other pyramid 481 static void pyramid_scale_gradient(pyramid_t* pyramid, pyramid_t* pC) 482 { 483 while (pyramid != NULL) 484 { 485 const int size = pyramid->rows * pyramid->cols; 486 scale_gradient(size, pyramid->Gx, pC->Gx); 487 scale_gradient(size, pyramid->Gy, pC->Gy); 488 pyramid = pyramid->next; 489 pC = pC->next; 490 } 491 } 492 493 494 // free memory allocated for the pyramid 495 static void pyramid_free(pyramid_t* pyramid) 496 { 497 while (pyramid) 498 { 499 if(pyramid->Gx != NULL) 500 { 501 free(pyramid->Gx); 502 pyramid->Gx = NULL; 503 } 504 if(pyramid->Gy != NULL) 505 { 506 free(pyramid->Gy); 507 pyramid->Gy = NULL; 508 } 509 pyramid_t* const next = pyramid->next; 510 pyramid->prev = NULL; 511 pyramid->next = NULL; 512 free(pyramid); 513 pyramid = next; 514 } 515 } 516 517 518 // allocate memory for the pyramid 519 static pyramid_t * pyramid_allocate(int cols, int rows) 520 { 521 pyramid_t* level = NULL; 522 pyramid_t* pyramid = NULL; 523 pyramid_t* prev = NULL; 524 525 while(rows >= PYRAMID_MIN_PIXELS && cols >= PYRAMID_MIN_PIXELS) 526 { 527 level = (pyramid_t *) malloc(sizeof(pyramid_t)); 528 if(level == NULL) 529 { 530 fprintf(stderr, "ERROR: malloc in pyramid_alloc() (size:%d)", (int)sizeof(pyramid_t)); 531 exit(155); 532 } 533 memset( level, 0, sizeof(pyramid_t) ); 534 535 level->rows = rows; 536 level->cols = cols; 537 const int size = level->rows * level->cols; 538 level->Gx = matrix_alloc(size); 539 level->Gy = matrix_alloc(size); 540 541 level->prev = prev; 542 if(prev != NULL) 543 prev->next = level; 544 prev = level; 545 546 if(pyramid == NULL) 547 pyramid = level; 548 549 rows /= 2; 550 cols /= 2; 551 } 552 553 return pyramid; 554 } 555 556 557 // calculate gradients 558 static inline void calculate_gradient(const int cols, const int rows, const float* const lum, float* const Gx, float* const Gy) 559 { 560 #pragma omp parallel for schedule(static) 561 for(int ky=0; ky<rows; ky++){ 562 for(int kx=0; kx<cols; kx++){ 563 564 const int idx = kx + ky*cols; 565 566 if(kx == (cols - 1)) 567 Gx[idx] = 0; 568 else 569 Gx[idx] = lum[idx+1] - lum[idx]; 570 571 if(ky == (rows - 1)) 572 Gy[idx] = 0; 573 else 574 Gy[idx] = lum[idx + cols] - lum[idx]; 575 } 576 } 577 } 578 579 580 #define PFSEOL "\x0a" 581 static void dump_matrix_to_file(const int width, const int height, const float* const m, const char * const file_name) 582 { 583 FILE *fh = fopen( file_name, "wb" ); 584 // assert( fh != NULL ); 585 586 fprintf( fh, "PFS1" PFSEOL "%d %d" PFSEOL "1" PFSEOL "0" PFSEOL 587 "%s" PFSEOL "0" PFSEOL "ENDH", width, height, "Y"); 588 589 for( int y = 0; y < height; y++ ) 590 for( int x = 0; x < width; x++ ) { 591 int idx = x + y*width; 592 fwrite( &(m[idx]), sizeof( float ), 1, fh ); 593 } 594 595 fclose( fh ); 596 } 597 598 // calculate gradients for the pyramid 599 // lum_temp gets overwritten! 600 static void pyramid_calculate_gradient(pyramid_t* pyramid, float* lum_temp) 601 { 602 float* temp = matrix_alloc((pyramid->rows/2)*(pyramid->cols/2)); 603 float* const temp_saved = temp; 604 605 calculate_gradient(pyramid->cols, pyramid->rows, lum_temp, pyramid->Gx, pyramid->Gy); 606 607 pyramid = pyramid->next; 608 609 // int l = 1; 610 while(pyramid) 611 { 612 matrix_downsample(pyramid->prev->cols, pyramid->prev->rows, lum_temp, temp); 613 614 // fprintf( stderr, "%d x %d -> %d x %d\n", pyramid->cols, pyramid->rows, 615 // prev_pp->cols, prev_pp->rows ); 616 617 // char name[40]; 618 // sprintf( name, "ds_%d.pfs", l++ ); 619 // dump_matrix_to_file( pyramid->cols, pyramid->rows, temp, name ); 620 621 calculate_gradient(pyramid->cols, pyramid->rows, temp, pyramid->Gx, pyramid->Gy); 622 623 float* const dummy = lum_temp; 624 lum_temp = temp; 625 temp = dummy; 626 627 pyramid = pyramid->next; 628 } 629 630 matrix_free(temp_saved); 631 } 632 633 634 // x = -0.25 * b 635 static inline void solveX(const int n, const float* const b, float* const x) 636 { 637 #pragma omp parallel for schedule(static) 638 for (int i=0; i<n; i++) 639 x[i] = -0.25f * b[i]; 640 } 641 642 // divG_sum = A * x = sum(divG(x)) 643 // memory for the temporary pyramid px should be allocated 644 static inline void multiplyA(pyramid_t* px, pyramid_t* pC, const float* const x, float* const divG_sum) 645 { 646 matrix_copy(pC->rows*pC->cols, x, divG_sum); // use divG_sum as a temp variable 647 pyramid_calculate_gradient(px, divG_sum); 648 pyramid_scale_gradient(px, pC); // scale gradients by Cx,Cy from main pyramid 649 pyramid_calculate_divergence_sum(px, divG_sum); // calculate the sum of divergences 650 } 651 652 653 // bi-conjugate linear equation solver 654 // overwrites pyramid! 655 static void linbcg(pyramid_t* pyramid, pyramid_t* pC, float* const b, float* const x, const int itmax, const float tol, pfstmo_progress_callback progress_cb) 656 { 657 const int rows = pyramid->rows; 658 const int cols = pyramid->cols; 659 const int n = rows*cols; 660 const float tol2 = tol*tol; 661 662 float* const z = matrix_alloc(n); 663 float* const zz = matrix_alloc(n); 664 float* const p = matrix_alloc(n); 665 float* const pp = matrix_alloc(n); 666 float* const r = matrix_alloc(n); 667 float* const rr = matrix_alloc(n); 668 float* const x_save = matrix_alloc(n); 669 670 const float bnrm2 = matrix_DotProduct(n, b, b); 671 672 multiplyA(pyramid, pC, x, r); // r = A*x = divergence(x) 673 matrix_subtract(n, b, r); // r = b - r 674 float err2 = matrix_DotProduct(n, r, r); // err2 = r.r 675 676 // matrix_copy(n, r, rr); // rr = r 677 multiplyA(pyramid, pC, r, rr); // rr = A*r 678 679 float bkden = 0; 680 float saved_err2 = err2; 681 matrix_copy(n, x, x_save); 682 683 const float ierr2 = err2; 684 const float percent_sf = 100.0f/logf(tol2*bnrm2/ierr2); 685 686 int iter = 0; 687 bool reset = true; 688 int num_backwards = 0; 689 const int num_backwards_ceiling = 3; 690 for (; iter < itmax; iter++) 691 { 692 fprintf( stderr, "Iter %d\n", iter ); 693 if( progress_cb != NULL ) 694 progress_cb( (int) (logf(err2/ierr2)*percent_sf)); 695 696 solveX(n, r, z); // z = ~A(-1) * r = -0.25 * r 697 solveX(n, rr, zz); // zz = ~A(-1) * rr = -0.25 * rr 698 699 const float bknum = matrix_DotProduct(n, z, rr); 700 701 if(reset) 702 { 703 reset = false; 704 matrix_copy(n, z, p); 705 matrix_copy(n, zz, pp); 706 } 707 else 708 { 709 const float bk = bknum / bkden; // beta = ... 710 #pragma omp parallel for schedule(static) 711 for (int i = 0; i < n; i++) 712 { 713 p[i] = z[i] + bk * p[i]; 714 pp[i] = zz[i] + bk * pp[i]; 715 } 716 } 717 718 bkden = bknum; // numerato becomes the dominator for the next iteration 719 720 multiplyA(pyramid, pC, p, z); // z = A* p = divergence( p) 721 multiplyA(pyramid, pC, pp, zz); // zz = A*pp = divergence(pp) 722 723 const float ak = bknum / matrix_DotProduct(n, z, pp); // alfa = ... 724 #pragma omp parallel for schedule(static) 725 for(int i = 0 ; i < n ; i++ ) 726 { 727 r[i] -= ak * z[i]; // r = r - alfa * z 728 rr[i] -= ak * zz[i]; //rr = rr - alfa * zz 729 } 730 731 const float old_err2 = err2; 732 err2 = matrix_DotProduct(n, r, r); 733 734 // Have we gone unstable? 735 if (err2 > old_err2) 736 { 737 // Save where we've got to if it's the best yet 738 if (num_backwards == 0 && old_err2 < saved_err2) 739 { 740 saved_err2 = old_err2; 741 matrix_copy(n, x, x_save); 742 } 743 744 num_backwards++; 745 } 746 else 747 { 748 num_backwards = 0; 749 } 750 751 #pragma omp parallel for schedule(static) 752 for(int i = 0 ; i < n ; i++ ) 753 x[i] += ak * p[i]; // x = x + alfa * p 754 755 if (num_backwards > num_backwards_ceiling) 756 { 757 // Reset 758 reset = true; 759 num_backwards = 0; 760 761 // Recover saved value 762 matrix_copy(n, x_save, x); 763 764 // r = Ax 765 multiplyA(pyramid, pC, x, r); 766 767 // r = b - r 768 matrix_subtract(n, b, r); 769 770 // err2 = r.r 771 err2 = matrix_DotProduct(n, r, r); 772 saved_err2 = err2; 773 774 // rr = A*r 775 multiplyA(pyramid, pC, r, rr); 776 } 777 778 // fprintf(stderr, "iter:%d err:%f\n", iter+1, sqrtf(err2/bnrm2)); 779 if(err2/bnrm2 < tol2) 780 break; 781 } 782 783 // Use the best version we found 784 if (err2 > saved_err2) 785 { 786 err2 = saved_err2; 787 matrix_copy(n, x_save, x); 788 } 789 790 if (err2/bnrm2 > tol2) 791 { 792 // Not converged 793 if( progress_cb != NULL ) 794 progress_cb( (int) (logf(err2/ierr2)*percent_sf)); 795 if (iter == itmax) 796 fprintf(stderr, "\npfstmo_mantiuk06: Warning: Not converged (hit maximum iterations), error = %g (should be below %g).\n", sqrtf(err2/bnrm2), tol); 797 else 798 fprintf(stderr, "\npfstmo_mantiuk06: Warning: Not converged (going unstable), error = %g (should be below %g).\n", sqrtf(err2/bnrm2), tol); 799 } 800 else if (progress_cb != NULL) 801 progress_cb(100); 802 803 804 matrix_free(x_save); 805 matrix_free(p); 806 matrix_free(pp); 807 matrix_free(z); 808 matrix_free(zz); 809 matrix_free(r); 810 matrix_free(rr); 811 } 812 813 814 // conjugate linear equation solver 815 // overwrites pyramid! 816 static void lincg(pyramid_t* pyramid, pyramid_t* pC, const float* const b, float* const x, const int itmax, const float tol, pfstmo_progress_callback progress_cb) 817 { 818 const int rows = pyramid->rows; 819 const int cols = pyramid->cols; 820 const int n = rows*cols; 821 const float tol2 = tol*tol; 822 823 float* const x_save = matrix_alloc(n); 824 float* const r = matrix_alloc(n); 825 float* const p = matrix_alloc(n); 826 float* const Ap = matrix_alloc(n); 827 828 // bnrm2 = ||b|| 829 const float bnrm2 = matrix_DotProduct(n, b, b); 830 831 // r = b - Ax 832 multiplyA(pyramid, pC, x, r); 833 matrix_subtract(n, b, r); 834 float rdotr = matrix_DotProduct(n, r, r); // rdotr = r.r 835 836 // p = r 837 matrix_copy(n, r, p); 838 839 // Setup initial vector 840 float saved_rdotr = rdotr; 841 matrix_copy(n, x, x_save); 842 843 const float irdotr = rdotr; 844 const float percent_sf = 100.0f/logf(tol2*bnrm2/irdotr); 845 int iter = 0; 846 int num_backwards = 0; 847 const int num_backwards_ceiling = 3; 848 for (; iter < itmax; iter++) 849 { 850 if( progress_cb != NULL ) { 851 int ret = progress_cb( (int) (logf(rdotr/irdotr)*percent_sf)); 852 if( ret == PFSTMO_CB_ABORT && iter > 0 ) // User requested abort 853 break; 854 } 855 856 // Ap = A p 857 multiplyA(pyramid, pC, p, Ap); 858 859 // alpha = r.r / (p . Ap) 860 const float alpha = rdotr / matrix_DotProduct(n, p, Ap); 861 862 // r = r - alpha Ap 863 #pragma omp parallel for schedule(static) 864 for (int i = 0; i < n; i++) 865 r[i] -= alpha * Ap[i]; 866 867 // rdotr = r.r 868 const float old_rdotr = rdotr; 869 rdotr = matrix_DotProduct(n, r, r); 870 871 // Have we gone unstable? 872 if (rdotr > old_rdotr) 873 { 874 // Save where we've got to 875 if (num_backwards == 0 && old_rdotr < saved_rdotr) 876 { 877 saved_rdotr = old_rdotr; 878 matrix_copy(n, x, x_save); 879 } 880 881 num_backwards++; 882 } 883 else 884 { 885 num_backwards = 0; 886 } 887 888 // x = x + alpha p 889 #pragma omp parallel for schedule(static) 890 for (int i = 0; i < n; i++) 891 x[i] += alpha * p[i]; 892 893 894 // Exit if we're done 895 // fprintf(stderr, "iter:%d err:%f\n", iter+1, sqrtf(rdotr/bnrm2)); 896 if(rdotr/bnrm2 < tol2) 897 break; 898 899 if (num_backwards > num_backwards_ceiling) 900 { 901 // Reset 902 num_backwards = 0; 903 matrix_copy(n, x_save, x); 904 905 // r = Ax 906 multiplyA(pyramid, pC, x, r); 907 908 // r = b - r 909 matrix_subtract(n, b, r); 910 911 // rdotr = r.r 912 rdotr = matrix_DotProduct(n, r, r); 913 saved_rdotr = rdotr; 914 915 // p = r 916 matrix_copy(n, r, p); 917 } 918 else 919 { 920 // p = r + beta p 921 const float beta = rdotr/old_rdotr; 922 #pragma omp parallel for schedule(static) 923 for (int i = 0; i < n; i++) 924 p[i] = r[i] + beta*p[i]; 925 } 926 } 927 928 // Use the best version we found 929 if (rdotr > saved_rdotr) 930 { 931 rdotr = saved_rdotr; 932 matrix_copy(n, x_save, x); 933 } 934 935 if (rdotr/bnrm2 > tol2) 936 { 937 // Not converged 938 if( progress_cb != NULL ) 939 progress_cb( (int) (logf(rdotr/irdotr)*percent_sf)); 940 if (iter == itmax) 941 fprintf(stderr, "\npfstmo_mantiuk06: Warning: Not converged (hit maximum iterations), error = %g (should be below %g).\n", sqrtf(rdotr/bnrm2), tol); 942 else 943 fprintf(stderr, "\npfstmo_mantiuk06: Warning: Not converged (going unstable), error = %g (should be below %g).\n", sqrtf(rdotr/bnrm2), tol); 944 } 945 else if (progress_cb != NULL) 946 progress_cb(100); 947 948 matrix_free(x_save); 949 matrix_free(p); 950 matrix_free(Ap); 951 matrix_free(r); 952 } 953 954 955 // in_tab and out_tab should contain inccreasing float values 956 static inline float lookup_table(const int n, const float* const in_tab, const float* const out_tab, const float val) 957 { 958 if(unlikely(val < in_tab[0])) 959 return out_tab[0]; 960 961 for (int j = 1; j < n; j++) 962 if(val < in_tab[j]) 963 { 964 const float dd = (val - in_tab[j-1]) / (in_tab[j] - in_tab[j-1]); 965 return out_tab[j-1] + (out_tab[j] - out_tab[j-1]) * dd; 966 } 967 968 return out_tab[n-1]; 969 } 970 971 972 // transform gradient (Gx,Gy) to R 973 static inline void transform_to_R(const int n, float* const G) 974 { 975 #pragma omp parallel for schedule(static) 976 for(int j=0;j<n;j++) 977 { 978 // G to W 979 const float absG = fabsf(G[j]); 980 int sign; 981 if(G[j] < 0) 982 sign = -1; 983 else 984 sign = 1; 985 G[j] = (powf(10,absG) - 1.0f) * sign; 986 987 // W to RESP 988 if(G[j] < 0) 989 sign = -1; 990 else 991 sign = 1; 992 993 G[j] = sign * lookup_table(LOOKUP_W_TO_R, W_table, R_table, fabsf(G[j])); 994 } 995 } 996 997 // transform gradient (Gx,Gy) to R for the whole pyramid 998 static inline void pyramid_transform_to_R(pyramid_t* pyramid) 999 { 1000 while (pyramid != NULL) 1001 { 1002 const int size = pyramid->rows * pyramid->cols; 1003 transform_to_R(size, pyramid->Gx); 1004 transform_to_R(size, pyramid->Gy); 1005 pyramid = pyramid->next; 1006 } 1007 } 1008 1009 // transform from R to G 1010 static inline void transform_to_G(const int n, float* const R){ 1011 1012 #pragma omp parallel for schedule(static) 1013 for(int j=0;j<n;j++){ 1014 1015 // RESP to W 1016 int sign; 1017 if(R[j] < 0) 1018 sign = -1; 1019 else 1020 sign = 1; 1021 R[j] = sign * lookup_table(LOOKUP_W_TO_R, R_table, W_table, fabsf(R[j])); 1022 1023 // W to G 1024 if(R[j] < 0) 1025 sign = -1; 1026 else 1027 sign = 1; 1028 R[j] = log10f(fabsf(R[j]) + 1.0f) * sign; 1029 1030 } 1031 } 1032 1033 // transform from R to G for the pyramid 1034 static inline void pyramid_transform_to_G(pyramid_t* pyramid) 1035 { 1036 while (pyramid != NULL) 1037 { 1038 transform_to_G(pyramid->rows*pyramid->cols, pyramid->Gx); 1039 transform_to_G(pyramid->rows*pyramid->cols, pyramid->Gy); 1040 pyramid = pyramid->next; 1041 } 1042 } 1043 1044 // multiply gradient (Gx,Gy) values by float scalar value for the whole pyramid 1045 static inline void pyramid_gradient_multiply(pyramid_t* pyramid, const float val) 1046 { 1047 while (pyramid != NULL) 1048 { 1049 matrix_multiply_const(pyramid->rows*pyramid->cols, pyramid->Gx, val); 1050 matrix_multiply_const(pyramid->rows*pyramid->cols, pyramid->Gy, val); 1051 pyramid = pyramid->next; 1052 } 1053 } 1054 1055 1056 static int sort_float(const void* const v1, const void* const v2) 1057 { 1058 if (*((float*)v1) < *((float*)v2)) 1059 return -1; 1060 1061 if(likely(*((float*)v1) > *((float*)v2))) 1062 return 1; 1063 1064 return 0; 1065 } 1066 1067 1068 // transform gradients to luminance 1069 static void transform_to_luminance(pyramid_t* pp, float* const x, pfstmo_progress_callback progress_cb, const bool bcg, const int itmax, const float tol) 1070 { 1071 pyramid_t* pC = pyramid_allocate(pp->cols, pp->rows); 1072 pyramid_calculate_scale_factor(pp, pC); // calculate (Cx,Cy) 1073 pyramid_scale_gradient(pp, pC); // scale small gradients by (Cx,Cy); 1074 1075 float* b = matrix_alloc(pp->cols * pp->rows); 1076 pyramid_calculate_divergence_sum(pp, b); // calculate the sum of divergences (equal to b) 1077 1078 // calculate luminances from gradients 1079 if (bcg) 1080 linbcg(pp, pC, b, x, itmax, tol, progress_cb); 1081 else 1082 lincg(pp, pC, b, x, itmax, tol, progress_cb); 1083 1084 matrix_free(b); 1085 pyramid_free(pC); 1086 } 1087 1088 1089 struct hist_data 1090 { 1091 float size; 1092 float cdf; 1093 int index; 1094 }; 1095 1096 static int hist_data_order(const void* const v1, const void* const v2) 1097 { 1098 if (((struct hist_data*) v1)->size < ((struct hist_data*) v2)->size) 1099 return -1; 1100 1101 if (((struct hist_data*) v1)->size > ((struct hist_data*) v2)->size) 1102 return 1; 1103 1104 return 0; 1105 } 1106 1107 1108 static int hist_data_index(const void* const v1, const void* const v2) 1109 { 1110 return ((struct hist_data*) v1)->index - ((struct hist_data*) v2)->index; 1111 } 1112 1113 1114 static void contrast_equalization( pyramid_t *pp, const float contrastFactor ) 1115 { 1116 // Count sizes 1117 int total_pixels = 0; 1118 pyramid_t* l = pp; 1119 while (l != NULL) 1120 { 1121 total_pixels += l->rows * l->cols; 1122 l = l->next; 1123 } 1124 1125 // Allocate memory 1126 struct hist_data* hist = (struct hist_data*) malloc(sizeof(struct hist_data) * total_pixels); 1127 if (hist == NULL) 1128 { 1129 fprintf(stderr, "ERROR: malloc in contrast_equalization() (size:%d)", (int)sizeof(struct hist_data) * total_pixels); 1130 exit(155); 1131 } 1132 1133 // Build histogram info 1134 l = pp; 1135 int index = 0; 1136 while( l != NULL ) 1137 { 1138 const int pixels = l->rows*l->cols; 1139 const int offset = index; 1140 #pragma omp parallel for schedule(static) 1141 for(int c = 0; c < pixels; c++) 1142 { 1143 hist[c+offset].size = sqrtf( l->Gx[c]*l->Gx[c] + l->Gy[c]*l->Gy[c] ); 1144 hist[c+offset].index = c + offset; 1145 } 1146 index += pixels; 1147 l = l->next; 1148 } 1149 1150 // Generate histogram 1151 qsort(hist, total_pixels, sizeof(struct hist_data), hist_data_order); 1152 1153 // Calculate cdf 1154 const float norm = 1.0f / (float) total_pixels; 1155 #pragma omp parallel for schedule(static) 1156 for (int i = 0; i < total_pixels; i++) 1157 hist[i].cdf = ((float) i) * norm; 1158 1159 // Recalculate in terms of indexes 1160 qsort(hist, total_pixels, sizeof(struct hist_data), hist_data_index); 1161 1162 //Remap gradient magnitudes 1163 l = pp; 1164 index = 0; 1165 while( l != NULL ) { 1166 const int pixels = l->rows*l->cols; 1167 const int offset = index; 1168 1169 #pragma omp parallel for schedule(static) 1170 for( int c = 0; c < pixels; c++) 1171 { 1172 const float scale = contrastFactor * hist[c+offset].cdf/hist[c+offset].size; 1173 l->Gx[c] *= scale; 1174 l->Gy[c] *= scale; 1175 } 1176 index += pixels; 1177 l = l->next; 1178 } 1179 1180 free(hist); 1181 } 1182 1183 1184 // tone mapping 1185 int tmo_mantiuk06_contmap(const int c, const int r, float* const R, float* const G, float* const B, float* const Y, const float contrastFactor, const float saturationFactor, const bool bcg, const int itmax, const float tol, pfstmo_progress_callback progress_cb) 1186 { 1187 1188 const int n = c*r; 1189 1190 /* Normalize */ 1191 float Ymax = Y[0]; 1192 for (int j = 1; j < n; j++) 1193 if (Y[j] > Ymax) 1194 Ymax = Y[j]; 1195 1196 const float clip_min = 1e-7f*Ymax; 1197 #pragma omp parallel for schedule(static) 1198 for (int j = 0; j < n; j++) 1199 { 1200 if( unlikely(R[j] < clip_min) ) R[j] = clip_min; 1201 if( unlikely(G[j] < clip_min) ) G[j] = clip_min; 1202 if( unlikely(B[j] < clip_min) ) B[j] = clip_min; 1203 if( unlikely(Y[j] < clip_min) ) Y[j] = clip_min; 1204 } 1205 1206 #pragma omp parallel for schedule(static) 1207 for(int j=0;j<n;j++) 1208 { 1209 R[j] /= Y[j]; 1210 G[j] /= Y[j]; 1211 B[j] /= Y[j]; 1212 Y[j] = log10f(Y[j]); 1213 } 1214 1215 pyramid_t* pp = pyramid_allocate(c,r); // create pyramid 1216 float* tY = matrix_alloc(n); 1217 matrix_copy(n, Y, tY); // copy Y to tY 1218 pyramid_calculate_gradient(pp,tY); // calculate gradients for pyramid, destroys tY 1219 matrix_free(tY); 1220 pyramid_transform_to_R(pp); // transform gradients to R 1221 1222 /* Contrast map */ 1223 if( contrastFactor > 0.0f ) 1224 pyramid_gradient_multiply(pp, contrastFactor); // Contrast mapping 1225 else 1226 contrast_equalization(pp, -contrastFactor); // Contrast equalization 1227 1228 pyramid_transform_to_G(pp); // transform R to gradients 1229 transform_to_luminance(pp, Y, progress_cb, bcg, itmax, tol); // transform gradients to luminance Y 1230 pyramid_free(pp); 1231 1232 /* Renormalize luminance */ 1233 float* temp = matrix_alloc(n); 1234 1235 matrix_copy(n, Y, temp); // copy Y to temp 1236 qsort(temp, n, sizeof(float), sort_float); // sort temp in ascending order 1237 1238 // const float median = (temp[(int)((n-1)/2)] + temp[(int)((n-1)/2+1)]) * 0.5f; // calculate median 1239 const float CUT_MARGIN = 0.1f; 1240 1241 float trim = (n-1) * CUT_MARGIN * 0.01f; 1242 float delta = trim - floorf(trim); 1243 const float l_min = temp[(int)floorf(trim)] * delta + temp[(int)ceilf(trim)] * (1.0f-delta); 1244 1245 trim = (n-1) * (100.0f - CUT_MARGIN) * 0.01f; 1246 delta = trim - floorf(trim); 1247 const float l_max = temp[(int)floorf(trim)] * delta + temp[(int)ceilf(trim)] * (1.0f-delta); 1248 1249 matrix_free(temp); 1250 1251 const float disp_dyn_range = 2.3f; 1252 #pragma omp parallel for schedule(static) 1253 for(int j=0;j<n;j++) 1254 Y[j] = (Y[j] - l_min) / (l_max - l_min) * disp_dyn_range - disp_dyn_range; // x scaled 1255 // 1256 /* Transform to linear scale RGB */ 1257 #pragma omp parallel for schedule(static) 1258 for(int j=0;j<n;j++) 1259 { 1260 Y[j] = powf(10,Y[j]); 1261 R[j] = powf( R[j], saturationFactor) * Y[j]; 1262 G[j] = powf( G[j], saturationFactor) * Y[j]; 1263 B[j] = powf( B[j], saturationFactor) * Y[j]; 1264 } 1265 1266 return PFSTMO_OK; 1267 } 1268