1 /**
2 * @file tmo_ferradans11.cpp
3 * Implementation of the algorithm presented in :
4 *
5 * An Analysis of Visual Adaptation and Contrast Perception for Tone Mapping
6 * S. Ferradans, M. Bertalmio, E. Provenzi, V. Caselles
7 * In IEEE Trans. Pattern Analysis and Machine Intelligence
8 *
9 *
10 * @author Sira Ferradans Copyright (C) 2013
11 *
12 *
13 * This file is a part of PFSTMO package.
14 * ----------------------------------------------------------------------
15 *
16 * This program is free software; you can redistribute it and/or modify
17 * it under the terms of the GNU General Public License as published by
18 * the Free Software Foundation; either version 2 of the License, or
19 * (at your option) any later version.
20 *
21 * This program is distributed in the hope that it will be useful,
22 * but WITHOUT ANY WARRANTY; without even the implied warranty of
23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 * GNU General Public License for more details.
25 *
26 * You should have received a copy of the GNU General Public License
27 * along with this program; if not, write to the Free Software
28 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
29 * ----------------------------------------------------------------------
30 *
31 */
32
33 #include <config.h>
34
35 #include <vector>
36 #include <algorithm>
37
38 #include <fftw3.h>
39 #include <assert.h>
40 #include <pfs.h>
41
42 #include "Imagen.h"
43 #include "pfstmo.h"
44 #include <pfs.h>
45
46 #include <cstring>
47 #include <iostream>
48
49 #include <stdlib.h>
50
51 using namespace std;
52
53
54
55 //for debugging purposes
56 #if 0
57 #define PFSEOL "\x0a"
58 static void dumpPFS( const char *fileName, const pfstmo::Array2D *data, const char *channelName )
59 {
60 FILE *fh = fopen( fileName, "wb" );
61 assert( fh != NULL );
62
63 int width = data->getCols();
64 int height = data->getRows();
65
66 fprintf( fh, "PFS1" PFSEOL "%d %d" PFSEOL "1" PFSEOL "0" PFSEOL
67 "%s" PFSEOL "0" PFSEOL "ENDH", width, height, channelName );
68
69 for( int y = 0; y < height; y++ )
70 for( int x = 0; x < width; x++ ) {
71 float d = (*data)(x,y);
72 fwrite( &d, sizeof( float ), 1, fh );
73 }
74
75 fclose( fh );
76 }
77 #endif
78
79
80 //--------------------------------------------------------------------
81
82
83 /*Implementation, hardcoded of the R function */
apply_arctg_slope10(float Ip,float I,float I2,float I3,float I4,float I5,float I6,float I7)84 float apply_arctg_slope10(float Ip,float I,float I2,float I3,float I4,float I5,float I6,float I7)
85 {
86 //Arctan, slope 10
87 float gr1,gr2,gr3,gr4,gr5,gr6,gr7;
88 gr1=Ip-I;
89 gr2=Ip*Ip-2*Ip*I+I2;
90 gr3=Ip*Ip*Ip-3*Ip*Ip*I+3*Ip*I2-I3;
91 gr4=Ip*Ip*Ip*Ip+4*Ip*Ip*I2+I4-4*Ip*I3+2*Ip*Ip*(I2-2*Ip*I);
92 gr5=Ip*gr4-Ip*Ip*Ip*Ip*I-4*Ip*Ip*I3-I5+4*Ip*I4-2*Ip*Ip*(I3-2*Ip*I2);
93 gr6=Ip*Ip*gr4-Ip*Ip*Ip*Ip*Ip*I-4*Ip*Ip*Ip*I3-Ip*I5+4*Ip*Ip*I4-2*Ip*Ip*Ip*(I3-2*Ip*I2)-(I*Ip*gr4-Ip*Ip*Ip*Ip*I2-4*Ip*Ip*I4-I6+4*Ip*I5-2*Ip*Ip*(I4-2*Ip*I3) );
94 gr7=Ip*Ip*(Ip*( Ip*Ip*Ip*Ip+4*Ip*Ip*I2+I4-4*Ip*I3+2*Ip*Ip*(I2-2*Ip*I) )-Ip*Ip*Ip*Ip*I-4*Ip*Ip*I3-I5+4*Ip*I4-2*Ip*Ip*(I3-2*Ip*I2))-2*Ip*(Ip*(Ip*Ip*Ip*Ip*I+4*Ip*Ip*I3+I5-4*Ip*I4+2*Ip*Ip*(I3-2*Ip*I2)) -Ip*Ip*Ip*Ip*I2-4*Ip*Ip*I4-I6+4*Ip*I5-2*Ip*Ip*(I4-2*Ip*I3))+Ip*(Ip*Ip*Ip*Ip*I2+4*Ip*Ip*I4+I6-4*Ip*I5+2*Ip*Ip*(I4-2*Ip*I3)) -Ip*Ip*Ip*Ip*I3-4*Ip*Ip*I5-I7+4*Ip*I6-2*Ip*Ip*(I5-2*Ip*I4);
95
96 // Hardcoded coefficients of the polynomial of degree 9 that allows to approximate the neighborhood averaging as a sum of convolutions
97 return(( -7.7456e+00)*gr7+ (3.1255e-16)*gr6+(1.5836e+01)*gr5+(-1.8371e-15)*gr4+(-1.1013e+01)*gr3+(4.4531e-16)*gr2+(3.7891e+00)*gr1+ 1.2391e-15 ) ;
98 }
99
100
101
102 /*
103 * This Quickselect routine is based on the algorithm described in
104 * "Numerical recipes in C", Second Edition,
105 * Cambridge University Press, 1992, Section 8.5, ISBN 0-521-43108-5
106 * This code by Nicolas Devillard - 1998. Public domain.
107 */
108
109
110 #define ELEM_SWAP(a,b) { register float t=(a);(a)=(b);(b)=t; }
111 //#define ELEM_SWAP(a,b) { register float t=a;a=b;b=t; }
112
quick_select(float arr[],int n)113 double quick_select(float arr[], int n)
114 {
115 int low, high ;
116 int median;
117 int middle, ll, hh;
118
119 low = 0 ; high = n-1 ; median = (low + high) / 2;
120 for (;;) {
121 if (high <= low) /* One element only */
122 return arr[median] ;
123
124 if (high == low + 1) { /* Two elements only */
125 if (arr[low] > arr[high])
126 ELEM_SWAP(arr[low], arr[high]) ;
127 return arr[median] ;
128 }
129
130 /* Find median of low, middle and high items; swap into position low */
131 middle = (low + high) / 2;
132 if (arr[middle] > arr[high]) ELEM_SWAP(arr[middle], arr[high]) ;
133 if (arr[low] > arr[high]) ELEM_SWAP(arr[low], arr[high]) ;
134 if (arr[middle] > arr[low]) ELEM_SWAP(arr[middle], arr[low]) ;
135
136 /* Swap low item (now in position middle) into position (low+1) */
137 ELEM_SWAP(arr[middle], arr[low+1]) ;
138
139 /* Nibble from each end towards middle, swapping items when stuck */
140 ll = low + 1;
141 hh = high;
142 for (;;) {
143 do ll++; while (arr[low] > arr[ll]) ;
144 do hh--; while (arr[hh] > arr[low]) ;
145
146 if (hh < ll)
147 break;
148
149 ELEM_SWAP(arr[ll], arr[hh]) ;
150 }
151
152 /* Swap middle item (in position low) back into correct position */
153 ELEM_SWAP(arr[low], arr[hh]) ;
154
155 /* Re-set active partition */
156 if (hh <= median)
157 low = ll;
158 if (hh >= median)
159 high = hh - 1;
160 }
161 }
162
163 #undef ELEM_SWAP
164
tmo_ferradans11(int col,int fil,float * imR,float * imG,float * imB,float rho,float invalpha)165 void tmo_ferradans11(int col,int fil,float* imR, float* imG, float* imB,float rho, float invalpha){
166
167 int length=fil*col;
168 int colors=3;
169 float lambda;
170 float dt=0.2;//1e-1;//
171 float threshold_diff=dt/20.0;//1e-5;//
172
173
174 Imagen RGBorig[3],*pIm;
175
176 pIm = new Imagen(fil,col);
177 for (int c =0;c<3;c++){
178 RGBorig[c] = *pIm;
179
180 }
181
182 for(int i=0;i<length;i++){
183 RGBorig[0].datos[i] = (double)max(imR[i],0);
184 RGBorig[1].datos[i] = (double)max(imG[i],0);
185 RGBorig[2].datos[i] = (double)max(imB[i],0);
186 }
187
188 delete(pIm);
189
190 ///////////////////////////////////////////
191
192 char* aux_char = new char[1000];
193
194 Imagen RGB[3];
195
196 int iteration=0;
197 char indice[3];
198 float difference0=2000.0;
199 float difference=1000.0;
200
201 float med[3];
202 float valmax[3];
203 float valmin[3];
204
205 Imagen aux;
206 float median,mu[3];
207
208 int cInit=clock(); //start counting time for first step
209
210 fprintf(stderr,"inv_alpha=%f, rho=%f\n",invalpha,rho);
211
212 for(int k=0;k<3;k++)
213 {
214 RGBorig[k]+=1e-6;
215 aux=RGBorig[k];
216 median=quick_select(aux.datos,length);
217 mu[k]=pow(aux.medval(),0.5)*pow(median,0.5);
218
219 RGB[k]=RGBorig[k] ;
220 }
221
222 // SEMISATURATION CONSTANT SIGMA DEPENDS ON THE ILUMINATION
223 // OF THE BACKGROUND. DATA IN TABLA1 FROM VALETON+VAN NORREN.
224 // TAKE AS REFERENCE THE CHANNEL WITH MAXIMUM ILUMINATION.
225 float muMax=max(max(mu[0],mu[1]),mu[2]);
226 for(int k=0;k<3;k++)
227 {
228 //MOVE log(mu) EQUALLY FOR THE 3 COLOR CHANNELS: rho DOES NOT CHANGE
229 //PARAMETER OF OUR ALGORITHM
230 float x=log10(muMax)-log10(mu[k]);
231
232 float z= - 0.37*(log10(mu[k])+4-rho) + 1.9;
233 mu[k]*=pow(10,z);
234 }
235
236
237 // VALETON + VAN NORREN:
238 // range = 4 orders; r is half the range
239 // n=0.74 : EXPONENT in NAKA-RUSHTON formula
240 float r=2;
241 float n=0.74;
242 for(int k=0;k<3;k++)
243 {
244 //find ctes. for WEBER-FECHNER from NAKA-RUSHTON
245 float logs=log10(mu[k]);
246 float I0=mu[k]/pow(10,1.2);
247 float sigma_n=pow(mu[k],n);
248
249 // WYSZECKI-STILES, PÁG. 530: FECHNER FRACTION
250 float K_=100.0/1.85;
251 if(k==2)
252 K_= 100.0/8.7;
253 float Ir=pow(10,logs+r);
254 float mKlogc=pow(Ir,n)/(pow(Ir,n)+pow(mu[k],n))-K_*log10(Ir+I0);
255
256 //mix W-F and N-R
257 for(int d=0;d<length;d++)
258 {
259 float x=log10(RGBorig[k].datos[d]);
260 float In= pow(RGBorig[k].datos[d],n);
261 float srn=pow(pow(10,logs+r),n);
262 // before logs+r apply W-F, after N-R
263 if(x<=logs+r)
264 RGB[k].datos[d]=K_*log10( RGBorig[k].datos[d] + I0)+mKlogc;
265 else
266 RGB[k].datos[d]= In/(In+sigma_n);
267
268 }
269
270 float minmez=RGB[k].minval();
271 RGB[k]+= -minmez;
272 float escalamez=1.0/(RGB[k].maxval()+1e-12);
273 RGB[k]*=escalamez;
274
275 }
276
277 int cFin=clock();
278 fprintf(stderr,"Step 1 done in: %f", (float)(cFin - cInit)/(float)CLOCKS_PER_SEC);
279
280 for(int color=0;color<colors;color++)
281 {
282 RGBorig[color]=RGB[color];
283 med[color]=RGB[color].medval();
284 }
285
286 lambda=1;
287
288 int c1,c0=clock();
289
290 Imagen RGB0=RGB[0];
291 Imagen u=RGB[0];
292 Imagen u0=RGB[0];
293 Imagen u2=RGB[0];
294 Imagen u3=RGB[0];
295 Imagen u4=RGB[0];
296 Imagen u5=RGB[0];
297 Imagen u6=RGB[0];
298 Imagen u7=RGB[0];
299 Imagen I0= RGB0;
300 Imagen ut= RGB0;
301
302
303 fftwf_complex* U = (fftwf_complex*) fftw_malloc(sizeof(fftwf_complex) * length);
304 fftwf_plan pU = fftwf_plan_dft_r2c_2d(fil, col, u0.datos, U,FFTW_ESTIMATE);
305 fftwf_complex* U2 = (fftwf_complex*) fftw_malloc(sizeof(fftwf_complex) * length);
306 fftwf_plan pU2 = fftwf_plan_dft_r2c_2d(fil, col, u2.datos, U2,FFTW_ESTIMATE);
307 fftwf_complex* U3 = (fftwf_complex*) fftw_malloc(sizeof(fftwf_complex) * length);
308 fftwf_plan pU3 = fftwf_plan_dft_r2c_2d(fil, col, u3.datos, U3,FFTW_ESTIMATE);
309 fftwf_complex* U4 = (fftwf_complex*) fftw_malloc(sizeof(fftwf_complex) * length);
310 fftwf_plan pU4 = fftwf_plan_dft_r2c_2d(fil, col, u4.datos, U4,FFTW_ESTIMATE);
311
312 fftwf_complex* U5 = (fftwf_complex*) fftw_malloc(sizeof(fftwf_complex) * length);
313 fftwf_plan pU5 = fftwf_plan_dft_r2c_2d(fil, col, u5.datos, U5,FFTW_ESTIMATE);
314 fftwf_complex* U6 = (fftwf_complex*) fftw_malloc(sizeof(fftwf_complex) * length);
315 fftwf_plan pU6 = fftwf_plan_dft_r2c_2d(fil, col, u6.datos, U6,FFTW_ESTIMATE);
316 fftwf_complex* U7 = (fftwf_complex*) fftw_malloc(sizeof(fftwf_complex) * length);
317 fftwf_plan pU7 = fftwf_plan_dft_r2c_2d(fil, col, u7.datos, U7,FFTW_ESTIMATE);
318
319 fftwf_complex* UG = new fftwf_complex[fil*col];
320 fftwf_complex* U2G= new fftwf_complex[fil*col];
321 fftwf_complex* U3G= new fftwf_complex[fil*col];
322 fftwf_complex* U4G= new fftwf_complex[fil*col];
323 fftwf_complex* U5G= new fftwf_complex[fil*col];
324 fftwf_complex* U6G= new fftwf_complex[fil*col];
325 fftwf_complex* U7G= new fftwf_complex[fil*col];
326
327 Imagen iu(fil,col,0);
328 fftwf_plan pinvU = fftwf_plan_dft_c2r_2d(fil, col, UG,iu.datos, FFTW_ESTIMATE);
329 Imagen iu2(fil,col,0);
330 fftwf_plan pinvU2 = fftwf_plan_dft_c2r_2d(fil, col, U2G,iu2.datos, FFTW_ESTIMATE);
331 Imagen iu3(fil,col,0);
332 fftwf_plan pinvU3 = fftwf_plan_dft_c2r_2d(fil, col, U3G,iu3.datos, FFTW_ESTIMATE);
333 Imagen iu4(fil,col,0);
334 fftwf_plan pinvU4 = fftwf_plan_dft_c2r_2d(fil, col, U4G,iu4.datos, FFTW_ESTIMATE);
335 Imagen iu5(fil,col,0);
336 fftwf_plan pinvU5 = fftwf_plan_dft_c2r_2d(fil, col, U5G,iu5.datos, FFTW_ESTIMATE);
337 Imagen iu6(fil,col,0);
338 fftwf_plan pinvU6 = fftwf_plan_dft_c2r_2d(fil, col, U6G,iu6.datos, FFTW_ESTIMATE);
339 Imagen iu7(fil,col,0);
340 fftwf_plan pinvU7 = fftwf_plan_dft_c2r_2d(fil, col, U7G,iu7.datos, FFTW_ESTIMATE);
341
342 float alpha=min(fil,col)/invalpha;
343 Imagen g;
344 g=nucleo_gaussiano(fil,col,alpha);
345 g.escala(1.0,0.0);
346
347 g.fftshift();
348
349 float suma=0, norm=1.0/(fil*col);
350 for(int i=0;i<length;i++)
351 suma+=g.datos[i];
352
353 g*=(1.0/suma); //now integral of the weight sums 1
354
355 fftwf_complex* G = (fftwf_complex*) fftw_malloc(sizeof(fftwf_complex) * length);
356 fftwf_plan pG = fftwf_plan_dft_r2c_2d(fil, col, g.datos, G,FFTW_ESTIMATE);
357 fftwf_execute(pG);
358
359 while(difference>threshold_diff)
360 {
361 iteration++;
362 difference0=difference;
363 difference=0.0;
364 for(int color=0;color<colors;color++)
365 {
366 u0=RGB[color];
367 RGB0=RGB[color];
368
369 u=u0;
370 u2=u; u2*=u;
371 u3=u2;u3*=u;
372 u4=u3;u4*=u;
373 u5=u4;u5*=u;
374 u6=u5;u6*=u;
375 u7=u6;u7*=u;
376
377 fftwf_execute(pU);
378 fftwf_execute(pU2);
379 fftwf_execute(pU3);
380 fftwf_execute(pU4);
381 fftwf_execute(pU5);
382 fftwf_execute(pU6);
383 fftwf_execute(pU7);
384
385 producto(U,G,UG,fil,col);
386 producto(U2,G,U2G,fil,col);
387 producto(U3,G,U3G,fil,col);
388 producto(U4,G,U4G,fil,col);
389 producto(U5,G,U5G,fil,col);
390 producto(U6,G,U6G,fil,col);
391 producto(U7,G,U7G,fil,col);
392
393 fftwf_execute(pinvU);
394 iu *= norm;
395
396 fftwf_execute(pinvU2);
397 iu2 *= norm;
398 fftwf_execute(pinvU3);
399 iu3 *= norm;
400 fftwf_execute(pinvU4);
401 iu4 *= norm;
402 fftwf_execute(pinvU5);
403 iu5 *= norm;
404 fftwf_execute(pinvU6);
405 iu6 *= norm;
406 fftwf_execute(pinvU7);
407 iu7 *= norm;
408
409 for(int i=0;i<length;i++)
410 {
411 // compute contrast component
412 u0.datos[i]=apply_arctg_slope10(u0.datos[i],iu.datos[i],iu2.datos[i],iu3.datos[i],iu4.datos[i],iu5.datos[i],iu6.datos[i],iu7.datos[i]);
413
414 //project onto the interval [-1,1]
415 u0.datos[i] = max( min( u0.datos[i],1) , -1 );
416 }
417
418 // normalizing R term to estandarize results
419 u0 *= (1.0/u0.maxabsval());
420
421 double norm = (1.0 + dt*(1.0+255.0/253.0));// assuming alpha=255/253,beta=1
422 for(int i=0;i<length;i++)
423 {
424 RGB[color].datos[i] = (RGB[color].datos[i] + dt*(RGBorig[color].datos[i] + 0.5*u0.datos[i]+255.0/253.0 *med[color])) / norm;
425 //project onto the interval [0,1]
426 RGB[color].datos[i] = max( min( RGB[color].datos[i],1) , 0 );
427 }
428
429 difference+=RGB0.MSE(RGB[color],1.0);
430
431 //for debugging purposes only:
432 //fprintf(stderr,"diff=%f, it=%d, range:(%f,%f)\n",difference,iteration, RGB[color].minval(),RGB[color].maxval());
433 }
434 c1=clock();
435 }
436 delete[] U;
437 delete[] U2;
438 delete[] U3;
439 delete[] U4;
440 delete[] U5;
441 delete[] U6;
442 delete[] U7;
443
444 delete[] UG;
445 delete[] U2G;
446 delete[] U3G;
447 delete[] U4G;
448 delete[] U5G;
449 delete[] U6G;
450 delete[] U7G;
451
452 fprintf(stderr,"\nComplete execution done in: %f secs\n", (float)(c1 - cInit)/(float)CLOCKS_PER_SEC);
453
454 for(int c=0;c<3;c++)
455 RGB[c].escala(1,0);
456
457 //range between (0,1)
458 for (int i =0;i<length;i++){
459 imR[i] = (float)RGB[0].datos[i];
460 imG[i] = (float)RGB[1].datos[i];
461 imB[i] = (float)RGB[2].datos[i];
462 }
463
464 delete[] G;
465 }
466
467
468
469