← Back to branch summary

~ubuntu-branches/ubuntu/natty/qtpfsgui/natty

« back to all changes in this revision

Viewing changes to src/TM_operators/mantiuk06/contrast_domain.cpp

  • Committer: Bazaar Package Importer
  • Author(s): Cyril Brulebois
  • Date: 2008-01-06 04:39:36 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20080106043936-a9u9g7yih3w16ru5
Tags: 1.9.0-1
* New upstream release.
* Replace “COPYING” with “LICENSE” in the NOT_NEEDED variable of
  debian/rules, following upstream's renaming.
* Update debian/links accordingly.
* Delete the matching TODO item since there's no longer needed to have a
  patched (with HTML tags) license file to get a correct display in the
  “License agreement” tab.
* Update the gcc4.3 patch (drop the hunk touching src/Libpfs/pfs.cpp):
   - 20_gcc4.3_includes.
* Add a link from /usr/share/qtpfsgui/html to the HTML documentation
  under /usr/share/doc/qtpfsgui/html since the former is used at runtime
  to display the manual.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/TM_operators/mantiuk06/contrast_domain.cpp
29
29
 * 
30
30
 * @author Radoslaw Mantiuk, <radoslaw.mantiuk@gmail.com>
31
31
 *
32
 
 * $Id: contrast_domain.cpp,v 1.5 2007/06/16 19:23:08 rafm Exp $
 
32
 * $Id: contrast_domain.cpp,v 1.6 2007/10/17 10:03:21 rafm Exp $
33
33
 */
34
34
 
35
35
 
38
38
#include <string.h>
39
39
#include <math.h>
40
40
 
 
41
//#include <sys/time.h>
 
42
 
41
43
#include "contrast_domain.h"
42
44
 
43
45
typedef struct {
117
119
  float matrix_median(int n, float* m);
118
120
  void pyramid_gradient_multiply(PYRAMID* pyramid, float val);
119
121
 
 
122
  void dump_matrix_to_file(int width, int height, float* m, const char *file_name);
120
123
  void matrix_show(char* text, int rows, int cols, float* data);
121
124
  void pyramid_show(PYRAMID* pyramid);
122
125
 
193
196
// upsampled matrix is twice bigger in each direction than data[]
194
197
// res should be a pointer to allocated memory for bigger matrix
195
198
// nearest neighborhood method - should be changed!
 
199
// cols and rows are the dimmensions of the output matrix
196
200
void ContrastDomain::matrix_upsample(int cols, int rows, float* in, float* out){
197
201
 
198
 
  float dx = (float)cols / (float)(cols*2);
199
 
  float dy = (float)rows / (float)(rows*2);
200
 
 
201
 
  const int outRows = rows * 2;
202
 
  const int outCols = cols * 2;
203
 
 
204
 
  const float inRows = rows;
205
 
  const float inCols = cols;
206
 
 
 
202
  const int inRows = (rows)/2;
 
203
  const int inCols = (cols)/2;
 
204
 
 
205
  const int outRows = rows;
 
206
  const int outCols = cols;
 
207
 
 
208
  const float dx = (float)inCols / ((float)outCols);
 
209
  const float dy = (float)inRows / ((float)outRows);
 
210
  
207
211
  const float filterSize = 1;
208
212
  
209
213
  float sx, sy;
233
237
  return;
234
238
}
235
239
 
 
240
 
 
241
// void ContrastDomain::matrix_upsample(int cols, int rows, float* in, float* out)
 
242
// {
 
243
//   const int inRows = rows/2;
 
244
//   const int inCols = cols/2;
 
245
 
 
246
//   const int outRows = rows;
 
247
//   const int outCols = cols;
 
248
  
 
249
//   int sx, sy;
 
250
//   int x, y;
 
251
//   for( y = 0, sy = 0; sy < inRows; y+=2, sy++ ) {
 
252
//     float *inLine = in + sy*inCols;
 
253
//     float *outLine = out + y*outCols;
 
254
//     for( x = 0, sx = 0; sx < inCols; x += 2, sx++ ) {
 
255
//       float &v = inLine[sx];
 
256
//       outLine[x] = v;
 
257
//       outLine[x+1] = v;
 
258
//       outLine[x+outCols] = v;
 
259
//       outLine[x+outCols+1] = v;
 
260
//     }
 
261
//   }
 
262
//   // odd output size
 
263
//   if( outCols & 1 ) 
 
264
//     for( y = 0; y < (outRows& ~1); y++ )
 
265
//       out[y*outCols+outCols-1] = out[y*outCols+outCols-2];
 
266
  
 
267
//   if( outRows & 1 ) 
 
268
//     for( x = 0; x < outCols; x++ ) {
 
269
//       out[outCols*(outCols-1)+x] = out[outCols*(outCols-1)+x];
 
270
  
 
271
 
 
272
//   return;
 
273
// }
 
274
 
 
275
 
236
276
// downsample the matrix
237
277
void ContrastDomain::matrix_downsample(int cols, int rows, float* data, float* res){
238
278
 
239
 
  const float inRows = rows;
240
 
  const float inCols = cols;
241
 
 
242
 
  const int outRows = rows / 2;
243
 
  const int outCols = cols / 2;
244
 
 
245
 
  const float dx = (float)inCols / (float)outCols;
246
 
  const float dy = (float)inRows / (float)outRows;
 
279
      //=================
 
280
  const int inRows = rows;
 
281
  const int inCols = cols;
 
282
 
 
283
  const int outRows = (rows) / 2;
 
284
  const int outCols = (cols) / 2;
 
285
 
 
286
  const float dx = (float)inCols / ((float)outCols);
 
287
  const float dy = (float)inRows / ((float)outRows);
247
288
 
248
289
  const float filterSize = 0.5;
249
290
  
266
307
  return;
267
308
}
268
309
 
 
310
// void ContrastDomain::matrix_downsample(int cols, int rows, float* data, float* res){
 
311
 
 
312
//   const float inRows = rows;
 
313
//   const float inCols = cols;
 
314
 
 
315
//   const int outRows = rows / 2;
 
316
//   const int outCols = cols / 2;
 
317
 
 
318
//   const float dx = (float)inCols / (float)outCols;
 
319
//   const float dy = (float)inRows / (float)outRows;
 
320
 
 
321
//   const float filterSize = 0.5;
 
322
  
 
323
//   float sx, sy;
 
324
//   int x, y;
 
325
  
 
326
//   for( y = 0, sy = dy/2-0.5; y < outRows; y++, sy += dy )
 
327
//     for( x = 0, sx = dx/2-0.5; x < outCols; x++, sx += dx ) {
 
328
 
 
329
//       float pixVal = 0;
 
330
//       float w = 0;
 
331
//       for( float ix = max( 0, ceilf( sx-dx*filterSize ) ); ix <= min( floorf( sx+dx*filterSize ), inCols-1 ); ix++ )
 
332
//         for( float iy = max( 0, ceilf( sy-dx*filterSize ) ); iy <= min( floorf( sy+dx*filterSize), inRows-1 ); iy++ ) {
 
333
//           pixVal += data[(int)ix + (int)iy * (int)inCols];
 
334
//           w += 1;
 
335
//         }     
 
336
//       res[x + y * outCols] = pixVal/w;      
 
337
//     }        
 
338
        
 
339
//   return;
 
340
// }
 
341
 
 
342
 
269
343
// return = a + b
270
344
float* ContrastDomain::matrix_add(int n, float* a, float* b){
271
345
  for(int i=0; i<n; i++){
384
458
// calculate sum of divergences in the pyramid
385
459
//      divergences for a particular levels should be calculated earlier 
386
460
//      and set in PYRAMID_LEVEL::divG
387
 
void ContrastDomain::pyramid_calculate_divergence_sum_in(PYRAMID* pyramid, float* divG_sum){
 
461
void ContrastDomain::pyramid_calculate_divergence_sum(PYRAMID* pyramid, float* divG_sum){
388
462
        
389
463
  if(pyramid->next != NULL)     
390
 
    pyramid_calculate_divergence_sum_in((PYRAMID*)pyramid->next, divG_sum);
 
464
    pyramid_calculate_divergence_sum((PYRAMID*)pyramid->next, divG_sum);
391
465
 
392
466
  PYRAMID_LEVEL* pl = pyramid->level;
393
 
        
394
 
  matrix_add(pl->rows * pl->cols, pl->divG, divG_sum);
395
 
        
 
467
  
396
468
  float* temp = matrix_alloc(pl->rows*pl->cols);
397
 
  matrix_copy(pl->rows*pl->cols, divG_sum, temp);
 
469
  if(pyramid->next != NULL)
 
470
    matrix_upsample(pl->cols, pl->rows, divG_sum, temp );
 
471
  else
 
472
    matrix_zero(pl->rows * pl->cols, temp);
 
473
  
 
474
  matrix_add(pl->rows * pl->cols, pl->divG, temp);
398
475
 
399
 
  matrix_upsample(pl->cols, pl->rows, temp, divG_sum);
 
476
//   char name[256];
 
477
//   sprintf( name, "Up_%d.pfs", pl->cols );
 
478
//   dump_matrix_to_file( pl->cols, pl->rows, temp, name );  
 
479
  
 
480
  matrix_copy(pl->rows*pl->cols, temp, divG_sum);
400
481
        
401
482
  matrix_free(temp);
402
483
 
406
487
 
407
488
// calculate the sum of divergences for the all pyramid level
408
489
// the smaller divergence map is upsamled and added to the divergence map for the higher level of pyramid
409
 
void ContrastDomain::pyramid_calculate_divergence_sum(PYRAMID* pyramid, float* divG_sum){
410
 
        
411
 
  PYRAMID_LEVEL* pl = pyramid->level;
412
 
  matrix_zero(pl->rows * pl->cols, divG_sum);
413
 
        
414
 
  if(pyramid->next != NULL)     
415
 
    pyramid_calculate_divergence_sum_in((PYRAMID*)pyramid->next, divG_sum);
 
490
// void ContrastDomain::pyramid_calculate_divergence_sum(PYRAMID* pyramid, float* divG_sum){
 
491
        
 
492
//   PYRAMID_LEVEL* pl = pyramid->level;
 
493
// //  matrix_zero(pl->rows * pl->cols, divG_sum);
 
494
        
 
495
// //  if(pyramid->next != NULL)        
 
496
//   pyramid_calculate_divergence_sum_in((PYRAMID*)pyramid, divG_sum);
416
497
 
417
 
  matrix_add(pl->rows * pl->cols, pl->divG, divG_sum);
 
498
// //  matrix_add(pl->rows * pl->cols, pl->divG, divG_sum);
418
499
                
419
 
  return;
420
 
}
 
500
//   return;
 
501
// }
421
502
 
422
503
// calculate scale factors (Cx,Cy) for gradients (Gx,Gy)
423
504
// C is equal to EDGE_WEIGHT for gradients smaller than GFIXATE or 1.0 otherwise
429
510
  float* C = matrix_alloc(n);
430
511
 
431
512
  for(int i=0; i<n; i++){
432
 
 
433
 
    if(fabs(G[i]) < GFIXATE)
434
 
      C[i] = EDGE_WEIGHT;
435
 
    else        
436
 
      C[i] = 1.0;
 
513
    
 
514
//    const float detectT = 0.001;
 
515
//    const float g = max( detectT, fabs(G[i]) );    
 
516
//    const float a = 0.038737;
 
517
//    const float b = 0.537756;
 
518
//    C[i] = a*pow(g,b);
 
519
    
 
520
     if(fabs(G[i]) < GFIXATE)
 
521
       C[i] = EDGE_WEIGHT;
 
522
     else       
 
523
       C[i] = 1.0;
437
524
  }
438
525
 
439
526
  return C;
553
640
    if(p == NULL)
554
641
      p = pyramid;
555
642
                
556
 
    rows /= 2;
557
 
    cols /= 2;          
 
643
    rows = (rows)/2;
 
644
    cols = (cols)/2;            
558
645
    if(rows < PYRAMID_MIN_PIXELS || cols < PYRAMID_MIN_PIXELS)
559
646
      break;
560
647
  }
588
675
  return;
589
676
}
590
677
 
 
678
#define PFSEOL "\x0a"
 
679
void ContrastDomain::dump_matrix_to_file(int width, int height, float* m, const char *file_name){
 
680
 
 
681
  FILE *fh = fopen( file_name, "wb" );
 
682
  assert( fh != NULL );
 
683
  
 
684
  fprintf( fh, "PFS1" PFSEOL "%d %d" PFSEOL "1" PFSEOL "0" PFSEOL
 
685
    "%s" PFSEOL "0" PFSEOL "ENDH", width, height, "Y");
 
686
 
 
687
  for( int y = 0; y < height; y++ )
 
688
    for( int x = 0; x < width; x++ ) {
 
689
      int idx = x + y*width;
 
690
      fwrite( &(m[idx]), sizeof( float ), 1, fh );
 
691
    }
 
692
  
 
693
  fclose( fh );
 
694
}  
591
695
 
592
696
// calculate gradients for the pyramid
593
697
void ContrastDomain::pyramid_calculate_gradient(PYRAMID* pyramid, float* lum){
600
704
  calculate_gradient(pl->cols, pl->rows, lum_temp, pl->Gx, pl->Gy);     
601
705
 
602
706
  PYRAMID* pp = (PYRAMID*)pyramid->next;
 
707
  PYRAMID* prev_pp = pyramid;
603
708
  float* temp;
604
709
 
605
 
  while(1){
606
 
    if(pp == NULL)
607
 
      break;
 
710
//   int l = 1;
 
711
  while( pp != NULL ){
 
712
    
608
713
    pl = pp->level;
609
714
                        
610
715
    temp = matrix_alloc(pl->rows * pl->cols);                   
611
 
    matrix_downsample(pl->cols*2, pl->rows*2, lum_temp, temp);  
 
716
    matrix_downsample(prev_pp->level->cols, prev_pp->level->rows, lum_temp, temp);
 
717
 
 
718
//     fprintf( stderr, "%d x %d -> %d x %d\n", pl->cols, pl->rows,
 
719
//       prev_pp->level->cols, prev_pp->level->rows );
 
720
    
 
721
//      char name[40];
 
722
//      sprintf( name, "ds_%d.pfs", l++ );
 
723
//      dump_matrix_to_file( pl->cols, pl->rows, temp, name );    
612
724
                
613
725
    calculate_gradient(pl->cols, pl->rows, temp, pl->Gx, pl->Gy);
614
726
                
615
727
    matrix_free(lum_temp);
616
728
    lum_temp = temp;    
617
 
                        
 
729
 
 
730
    prev_pp = pp;
618
731
    pp = (PYRAMID*)pp->next;    
619
732
  }
620
733
  matrix_free(lum_temp);
627
740
void ContrastDomain::solveX(int n, float* b, float* x){
628
741
 
629
742
  matrix_copy(n, b, x); // x = b
630
 
  matrix_multiply_const(n, x, -0.25);
 
743
  matrix_multiply_const(n, x, -0.25);  
631
744
}
632
745
 
633
746
// divG_sum = A * x = sum(divG(x))
665
778
        
666
779
  float* x = matrix_alloc(n); // allocate x matrix filled with zeros
667
780
  matrix_zero(n, x); // x = 0
668
 
        
 
781
  
669
782
  multiplyA(pyramid_tmp, pyramid, x, r); // r = A*x = divergence(x)
670
783
 
671
784
  matrix_substract(n, b, r); // r = b - r
672
785
 
673
786
  matrix_copy(n, r, rr); // rr = r
674
 
        
 
787
 
 
788
  multiplyA(pyramid_tmp, pyramid, r, rr); // rr = A*r 
 
789
  
675
790
  float bnrm;
676
791
        
677
792
  bnrm = sqrt(matrix_DotProduct(n, b, b));
679
794
  solveX(n, r, z); // z = ~A(-1) * r = -0.25 * r
680
795
        
681
796
#define TOL 0.01
682
 
#define ITMAX 40
 
797
#define ITMAX 30
683
798
  int iter = 0;
684
799
 
685
800
  float bknum=0;
720
835
    akden = matrix_DotProduct(n, z, pp);  // alfa denominator
721
836
                
722
837
    ak = bknum / akden; // alfa = ...
723
 
                
 
838
 
 
839
    // This may need transpossing - possible bug
724
840
    multiplyA(pyramid_tmp, pyramid, pp, zz); // zz = A*pp = divergence(pp)
725
841
 
726
842
    // values for the next iterration
766
882
// in_tab and out_tab should contain inscresing float values
767
883
float ContrastDomain::lookup_table(int n, float* in_tab, float* out_tab, float val){
768
884
 
769
 
  float ret;
 
885
  float ret = out_tab[0];
770
886
  float dd;
771
887
 
772
888
  if(val < in_tab[0])
911
1027
 
912
1028
// transform gradients to luminance
913
1029
float* ContrastDomain::transform_to_luminance(PYRAMID* pp){
914
 
 
 
1030
  
915
1031
  pyramid_calculate_scale_factor(pp); // calculate (Cx,Cy)
916
1032
  pyramid_scale_gradient(pp); // scale small gradients by (Cx,Cy);
917
1033
  pyramid_calculate_divergence(pp); // calculate divergence for the pyramid
 
1034
//  dump_matrix_to_file( pp->next->level->cols, pp->next->level->rows, pp->next->level->Gx, "Gx.pfs" );
918
1035
        
919
1036
  float* b = matrix_alloc(pp->level->cols * pp->level->rows);
920
1037
  pyramid_calculate_divergence_sum(pp, b); // calculate the sum od divergences (equal to b)
921
 
        
 
1038
 
 
1039
//  dump_matrix_to_file( pp->level->cols, pp->level->rows, b, "B.pfs" );
 
1040
 
 
1041
//   PYRAMID *pi = pp;
 
1042
//   int i = 1;
 
1043
//   while( pi != NULL ) {
 
1044
//     PYRAMID_LEVEL* pl = pi->level;
 
1045
//     char name[256];
 
1046
//     sprintf( name, "Cx_%d.pfs", i );
 
1047
//     dump_matrix_to_file( pl->cols, pl->rows, pl->Cx, name );
 
1048
//     sprintf( name, "Cy_%d.pfs", i );
 
1049
//     dump_matrix_to_file( pl->cols, pl->rows, pl->Cx, name );
 
1050
//     pi = pi->next;
 
1051
//     i++;
 
1052
//   }
 
1053
 
922
1054
  float* x = linbcg(pp, b); // calculate luminances from gradients
923
1055
        
924
1056
  matrix_free(b);
1008
1140
        
1009
1141
  PYRAMID* pp = pyramid_allocate(c,r); // create pyramid
1010
1142
  pyramid_calculate_gradient(pp,Y); // calculate gradiens for pyramid
1011
 
        
1012
 
  pyramid_transform_to_R(pp); // transform gradients to R
1013
1143
 
 
1144
//  dump_matrix_to_file( pp->next->level->cols, pp->next->level->rows, pp->next->level->Gx, "Gx2.pfs" );
 
1145
  
 
1146
  pyramid_transform_to_R(pp); // transform gradients to R  
1014
1147
 
1015
1148
  if( contrastFactor != 0 ) {
1016
1149
    // Contrast mapping
1021
1154
  }
1022
1155
        
1023
1156
  pyramid_transform_to_G(pp); // transform R to gradients
1024
 
        
 
1157
  
1025
1158
  float* x = transform_to_luminance(pp); // transform gradients to luminance Y
1026
1159
        
1027
 
//      for (int i=0;i<rows*cols;i++) {
1028
 
//              fprintf(stderr,"v=%f\n", (*Yo)(i));
1029
 
//      }
1030
1160
  float* temp = matrix_alloc(n);
1031
1161
        
1032
1162
  matrix_copy(n, x, temp); // copy x to temp
1033
1163
  qsort(temp, n, sizeof(float), sort_median); // sort temp in ascending order
1034
1164
        
1035
 
  float median = (temp[(int)((n-1)/2)] + temp[(int)((n-1)/2+1)]) / 2.0; // calculate median
 
1165
//   float median = (temp[(int)((n-1)/2)] + temp[(int)((n-1)/2+1)]) / 2.0; // calculate median
1036
1166
  // c and r should be even
1037
1167
  float CUT_MARGIN = 0.1;
1038
1168
        
1045
1175
  float p_max = temp[(int)floor(trim)] * delta + temp[(int)ceil(trim)] * (1-delta);     
1046
1176
        
1047
1177
  matrix_free(temp);
1048
 
        
1049
1178
 
1050
 
  float d;
1051
 
  if( (p_max - median) > (median - p_min) )
1052
 
    d = p_max - median;
1053
 
  else 
1054
 
    d = median - p_min;
1055
 
        
1056
 
  float l_max = median + d;
1057
 
  float l_min = median - d;
1058
 
        
1059
 
  for(int j=0;j<n;j++){
1060
 
        
1061
 
    x[j] = (x[j] - l_min) / (l_max - l_min) * 2.5 - 2.5; // x scaled to <-2.5,0> range
 
1179
  const float disp_dyn_range = 2.3;
 
1180
  for(int j=0;j<n;j++) {        
 
1181
    x[j] = (x[j] - p_min) / (p_max - p_min) * disp_dyn_range - disp_dyn_range; // x scaled to <-2.5,0> range
1062
1182
  }
1063
1183
        
1064
1184
  for(int j=0;j<n;j++){ 
1084
1204
  ContrastDomain contrast;      
1085
1205
  contrast.tone_mapping( cols, rows, R, G, B, Y, contrastFactor, saturationFactor, progress_report );  
1086
1206
}
1087
 
 
1088
 
 
1089