"Fossies" - the Fresh Open Source Software Archive

Member "pfstools-2.2.0/src/tmo/mantiuk06/contrast_domain.cpp" (12 Aug 2021, 38487 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 "contrast_domain.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.

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