1 /**
2  * @brief Display Adaptive TMO
3  *
4  * From:
5  * Rafal Mantiuk, Scott Daly, Louis Kerofsky.
6  * Display Adaptive Tone Mapping.
7  * To appear in: ACM Transactions on Graphics (Proc. of SIGGRAPH'08) 27 (3)
8  * http://www.mpi-inf.mpg.de/resources/hdr/datmo/
9  *
10  * This file is a part of LuminanceHDR package, based on pfstmo.
11  * ----------------------------------------------------------------------
12  *  This program is free software; you can redistribute it and/or modify
13  *  it under the terms of the GNU General Public License as published by
14  *  the Free Software Foundation; either version 2 of the License, or
15  *  (at your option) any later version.
16  *
17  *  This program is distributed in the hope that it will be useful,
18  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
19  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
20  *  GNU General Public License for more details.
21  *
22  *  You should have received a copy of the GNU General Public License
23  *  along with this program; if not, write to the Free Software
24  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
25  * ----------------------------------------------------------------------
26  *
27  * @author Rafal Mantiuk, <mantiuk@gmail.com>
28  *
29  * $Id: display_function.cpp,v 1.3 2008/06/16 18:42:58 rafm Exp $
30  */
31 #include "../../sleef.c"
32 #include "../../opthelper.h"
33 #define pow_F(a,b) (xexpf(b*xlogf(a)))
34 
35 #include <boost/math/constants/constants.hpp>
36 
37 #include <assert.h>
38 #include <math.h>
39 #include <stdlib.h>
40 #include <string.h>
41 
42 #include "Libpfs/pfs.h"
43 
44 #include "arch/string.h"
45 #include "display_function.h"
46 
47 // ========== GGBA Display Function ==============
48 
DisplayFunctionGGBA(float gamma,float L_max,float L_black,float E_amb,float screen_refl)49 DisplayFunctionGGBA::DisplayFunctionGGBA(float gamma, float L_max,
50                                          float L_black, float E_amb,
51                                          float screen_refl) {
52     init(gamma, L_max, L_black, E_amb, screen_refl);
53 }
54 
init(float gamma,float L_max,float L_black,float E_amb,float screen_refl)55 void DisplayFunctionGGBA::init(float gamma, float L_max, float L_black,
56                                float E_amb, float screen_refl) {
57     this->gamma = gamma;
58     this->L_max = L_max;
59     this->L_black = L_black;
60     this->E_amb = E_amb;
61     this->screen_refl = screen_refl;
62     this->L_offset =
63         L_black + screen_refl / boost::math::double_constants::pi * E_amb;
64 }
65 
DisplayFunctionGGBA(const char * predefined)66 DisplayFunctionGGBA::DisplayFunctionGGBA(const char *predefined) {
67     if (!strcasecmp(predefined, "lcd_office")) {
68         init(2.2f, 100, 0.8f, 400, 0.01f);
69     } else if (!strcasecmp(predefined, "lcd")) {
70         init(2.2f, 200, 0.8f, 60, 0.01f);
71     } else if (!strcasecmp(predefined, "lcd_bright")) {
72         init(2.6f, 500, 0.5f, 10, 0.01f);
73     } else if (!strcasecmp(predefined, "crt")) {
74         init(2.2f, 80, 1.0f, 60, 0.02f);
75     } else {
76         throw pfs::Exception(
77             "Unknown display type. Recognized types: "
78             "lcd_office, lcd, lcd_bright, "
79             "crt.");
80     }
81 }
82 
print(FILE * fh)83 void DisplayFunctionGGBA::print(FILE *fh) {
84     fprintf(fh, "Display function: gamma-gain-black-ambient\n");
85     fprintf(fh, "   gamma = %g\t L_max = %g\t L_black = %g\n", (double)gamma,
86             (double)L_max, (double)L_black);
87     fprintf(fh, "   E_amb = %g\t k     = %g\n", (double)E_amb,
88             (double)screen_refl);
89 }
90 
inv_display(float L)91 float DisplayFunctionGGBA::inv_display(float L) {
92     if (L < L_offset) L = L_offset;
93     if (L > (L_offset + L_max)) L = L_offset + L_max;
94     return pow_F((L - L_offset) / (L_max - L_black), 1 / gamma);
95 }
96 
display(float pix)97 float DisplayFunctionGGBA::display(float pix) {
98     assert(pix >= 0 && pix <= 1);
99     return pow(pix, gamma) * (L_max - L_black) + L_offset;
100 }
101 
102 #ifdef __SSE2__
103 
inv_display(vfloat L)104 vfloat DisplayFunctionGGBA::inv_display(vfloat L) {
105     const vfloat voffset = F2V(L_offset);
106     const vfloat vmax = F2V(L_max);
107     L = vmaxf(L, voffset);
108     L = vminf(L, voffset + vmax);
109     return pow_F((L - voffset) / (vmax - F2V(L_black)),
110                       F2V(1.0f / gamma));
111 }
112 
display(vfloat pix)113 vfloat DisplayFunctionGGBA::display(vfloat pix) {
114     return pow_F(pix, F2V(gamma)) *
115                (F2V(L_max) - F2V(L_black)) +
116            F2V(L_offset);
117 }
118 
119 #endif
120 
121 // ========== LUT Display Function ==============
122 
123 static const int max_lut_size = 4096;
124 
DisplayFunctionLUT(const char * file_name)125 DisplayFunctionLUT::DisplayFunctionLUT(const char *file_name)
126     : pix_lut(NULL), L_lut(NULL) {
127     FILE *fh = fopen(file_name, "r");
128     if (fh == NULL) throw pfs::Exception("Cannot open lookup-table file");
129 
130     L_lut = new float[max_lut_size];
131     pix_lut = new float[max_lut_size];
132 
133     const size_t max_line = 100;
134     char buf[max_line];
135 
136     int i = 0;
137     while (fgets(buf, max_line, fh) != NULL) {
138         float p_buf, L_buf;
139         if (sscanf(buf, "%f%*[ ,;]%f\n", &p_buf, &L_buf) != 2) continue;
140         if (p_buf < 0 || p_buf > 1) {
141             fclose(fh);
142             throw pfs::Exception(
143                 "Improper LUT: pixel values must be from 0 to 1");
144         }
145 
146         if (L_buf <= 0) {
147             fclose(fh);
148             throw pfs::Exception(
149                 "Improper LUT: luminance must be greater than 0");
150         }
151         L_lut[i] = log10(L_buf);
152         pix_lut[i] = p_buf;
153         i++;
154         if (i >= max_lut_size) {
155             fclose(fh);
156             throw pfs::Exception("LUT too large (more than 4096 entries)");
157         }
158     }
159     lut_size = i;
160     if (pix_lut[0] != 0 || pix_lut[lut_size - 1] != 1) {
161         fclose(fh);
162         throw pfs::Exception(
163             "The first and last LUT entries for pixel value should be 0 and 1");
164     }
165     fclose(fh);
166 }
167 
~DisplayFunctionLUT()168 DisplayFunctionLUT::~DisplayFunctionLUT() {
169     delete[] pix_lut;
170     delete[] L_lut;
171 }
172 
bin_search_interp(float x,const float * lut_x,const float * lut_y,const int lutSize)173 inline float bin_search_interp(float x, const float *lut_x, const float *lut_y,
174                                const int lutSize) {
175     if (x <= lut_x[0]) return lut_y[0];
176     if (x >= lut_x[lutSize - 1]) return lut_y[lutSize - 1];
177 
178     size_t l = 0, r = lutSize;
179     while (true) {
180         size_t m = (l + r) / 2;
181         if (m == l) break;
182         if (x < lut_x[m])
183             r = m;
184         else
185             l = m;
186     }
187     return lut_y[l] + (lut_y[l + 1] - lut_y[l]) * (x - lut_x[l]);
188 }
189 
inv_display(float L)190 float DisplayFunctionLUT::inv_display(float L) {
191     return bin_search_interp(log10(L), L_lut, pix_lut, lut_size);
192 }
193 
display(float pix)194 float DisplayFunctionLUT::display(float pix) {
195     return pow(10, bin_search_interp(pix, pix_lut, L_lut, lut_size));
196 }
197 
print(FILE * fh)198 void DisplayFunctionLUT::print(FILE *fh) {
199     fprintf(fh, "Display function: lookup-table\n");
200     fprintf(fh, "   L_min = %g \tL_max = %g\n", (double)pow(10, L_lut[0]),
201             (double)pow(10, L_lut[lut_size - 1]));
202 }
203 
204 // ========== Command line parsing ==============
205 /*
206 DisplayFunction *createDisplayFunctionFromArgs( int &argc, char* argv[] )
207 {
208   DisplayFunction *df = 0;
209 
210     for( int i=1 ; i<argc; i++ )
211     {
212       if( !strcmp( argv[i], "--display-function" ) || !strcmp( argv[i], "-d" )
213 ) {
214         if( i+1 >= argc )
215           throw pfs::Exception( "missing display function specification" );
216 
217         float gamma = 2.2f, L_max = 100.f, L_black = 1.f, k = 0.01, E_amb = 50;
218         bool GGBA_model = true;
219         char *token;
220         token = strtok( argv[i+1], ":" );
221         while( token != NULL ) {
222           if( !strncmp( token, "pd=", 3 )  ) {
223               df = new DisplayFunctionGGBA( token+3 );
224               GGBA_model = false;
225               break;
226           } else if( !strncmp( token, "lut=", 4 )  ) {
227               df = new DisplayFunctionLUT( token+4 );
228               GGBA_model = false;
229               break;
230           } else if( !strncmp( token, "g=", 2 ) ) {
231             gamma = strtod( token+2, NULL );
232           } else if( !strncmp( token, "l=", 2 ) ) {
233             L_max = strtod( token+2, NULL );
234           } else if( !strncmp( token, "b=", 2 ) ) {
235             L_black = strtod( token+2, NULL );
236           } else if( !strncmp( token, "k=", 2 ) ) {
237             k = strtod( token+2, NULL );
238           } else if( !strncmp( token, "a=", 2 ) ) {
239             E_amb = strtod( token+2, NULL );
240           } else {
241             throw pfs::Exception( "Bad display type specification" );
242           }
243           token = strtok( NULL, ":" );
244         }
245         if( GGBA_model )
246           df = new DisplayFunctionGGBA( gamma, L_max, L_black, E_amb, k );
247 
248         removeCommandLineArg( argc, argv, i, 2 );
249         break;
250       }
251     }
252     return df;
253 }
254 
255 static void removeCommandLineArg( int &argc, char* argv[],
256   int firstArgToRemove, int numArgsToRemove )
257 {
258   assert( firstArgToRemove+numArgsToRemove <= argc );
259   if( argc-firstArgToRemove-numArgsToRemove > 0 ) {
260     for( int i = firstArgToRemove; i < argc-numArgsToRemove; i++ )
261       argv[i] = argv[i+numArgsToRemove];
262   }
263 
264   argc -= numArgsToRemove;
265 }
266 
267 */
268